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SUMMARY 


PROBLEM 

Develop  a  computer  model  for  calculating  the  vertical  and  horizontal 
directionality  of  low-frequency  ambient  noise  for  an  ocean  environment. 


RESULTS 

A  FORTRAN  computer  model  has  been  developed.  The  Research 
Ambient  Noise  Directionality  Model  (RANDI)  has  given  results  in  good 
agreement  with  measured  data  for  the  Pacific,  Atlantic,  and  the  Mediterranean. 


RECOMMENDATIONS 

RANDI  noise  calculations  should  be  compared  with  ambient  noise 
measurements  for  other  ocean  areas,  seasons,  and  noise-source  distributions 
as  data,  sufficiently  documented  for  validation  purposes,  become  available. 
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13  ABSTRACT 

Naval  Ship  Systems  Command 

Washington,  D  C.  20360 

A  computer  model  for  calculating  the  vertical  and  horizontal  directionality  of 
ambient  noise  and  some  of  the  results  which  have  been  obtained  are  discussed.  The 
model  considers  three  sources  of  surface-generated  anisotropic  noise,  one  source  of 
surface-generated  isotropic  noise  and  two  sources  of  volumetric  isotropic  noise.  Noise 
levels  are  obtained  tor  various  propagation  conditions  and  for  different  sensor  depths. 
Results  agree  well  with  experimental  data,  and  it  is  shown  that  the  model  can  be  a 
useful  tool  in  the  planning  of  ambient  noise  measurement  experiments  and  in  the 
analysis  of  results.  1  __ 
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INTRODUCTION 


The  Research  Ambient  Noise  Directionality  Model  (RANDI)  is  an  operating  general- 
purpose  FORTRAN  Program  that  calculates  and  displays,  via  CALCOMP  plots,  the  vertical 
and  horizontal  directionality  of  low-frequency  (between  10  and  500  Hz)  ambient  noise  for 
an  ocean  environment.  This  first-phase  program  will  be  used  in  the  design  of  ambient  noise 
measurement  experiments  and  in  surveillance  systems  analysis  studies.  A  variable-dependent 
approach  based  on  system  and  environmental  parameter  variance  makes  this  model  inappro¬ 
priate  for  synoptic  forecasting.  It  will  be  improved  and  validated  as  time,  data,  and  state-of- 
the-art  permit. 

BACKGROUND 

The  ambient  noise  at  any  particular  point  in  the  ocean  depends  upon  the  amount  and 
location  of  noise-generating  sources  and  the  nature  of  the  local  acoustic  propagation  condi¬ 
tions,  i.e,,  how  the  noise  gets  from  the  generator  to  the  measuring  sensor. 

Sound  in  the  ocean  travels  in  curved  or  refracted  paths  and  will  arrive  at  a  hydrophone 
from  various  vertical  angles,  depending  on  the  depth  of  the  source  and  the  hydrophone  and 
on  the  separation  range.  The  path  an  individual  sound  ray  travels  can  include  surface  reflec¬ 
tions  and  bottom  bounces,  as  illustrated  in  the  ray  diagram  (Fig.  1 ).  The  bottom  bounce 
paths  will  generally  have  the  highest  sound  energy  losses,  while  the  paths  which  do  not  make 
contact  with  either  the  surface  or  the  bottom  will  have  the  smallest  transmission  losses.  The 
vertical  angle  between  one  of  these  latter  rays  and  the  horizontal  is  generally  less  than  1  5  deg 
at  any  range,  while  the  rays  making  contact  with  the  surface  or  bottom  can  have  vertical 
angles  approaching  90  deg.  This  effectively  divides  the  sound  rays  into  three  distinct  groups: 
(1 )  rays  which  come  in  contact  with  the  bottom  and  arrive  at  angles  from  approximately  15 
to  90  deg  down  from  the  horizontal;  (2)  the  near-horizontal  or  SOFAR  Channel  rays,  which 
touch  neither  surface  nor  bottom  and,  (3)  the  rays  arriving  at  angles  above  1  5  deg  and  which 
originated  or  were  reflected  from  near  the  surface. 

There  are  many  sources  of  noise,  the  most  important  being  winds,  waves,  shipping, 
thermal  agitation  and  biological  and  seismic  activity.  The  noise  due  to  winds,  waves,  and 
shipping  originates  near  the  surface  and  arrives  at  the  sensor  along  paths  which  travel  near 
or  reflect  from  the  surface.  In  the  case  of  distant  sources,  the  noise  travels  by  way  of  the 
SOFAR  Channel  and  w'ill  arrive,  at  a  sensor  located  within  the  channel,  from  nearly  hori¬ 
zontal  angles.  Noise  from  near  sources  (within  a  few  hundred  miles)  will  arrive  at  angles 
closer  to  the  vertical  than  does  the  SOFAR  Channel  noise.  Thermal  noise  comes  from 
throughout  the  medium  and  can  arrive  at  any  angle  or  from  any  direction.  Seismic  noise 
would  be  expected  to  arrive  along  those  paths  which  come  in  contact  with  radiating  or  re¬ 
flecting  surfaces,  including  the  bottom  and  distant  seamounts.  Biological  noise,  however, 
originates  along  the  surface,  the  bottom,  and  throughout  the  medium. 
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The  net  effect  of  propagating— to  the  sensor— the  noises  from  such  a  highly  stratified 
distribution  of  noise  generators  is  a  noise  field  having  vertical  directionality,  i.e.,  noise  level 
dependent  on  vertical  arrival  angle,  as  illustrated  in  Fig.  2. 

Various  ambient  noise  models  have  been  proposed  which  include  part  of  the  major 
noise  sources  and  invoke  certain  simplifying  assumptions  to  produce  a  medium  less  hostile 
to  mathematical  description  in  an  attempt  to  render  the  problem  solvable  and  yet  obtain, 
meaningful  results.  Few  investigators  have  succeeded,  and  then  only  in  highly  specialized 
cases.  There  is  documentation  on  several  of  these  ambient  noise  directionality  models.  In 
the  opinion  of  the  author,  none  is  adequate  for  use  in  the  design  of  detailed  directional 
noise  measurement  experiments  or  for  use  in  surveillance  systems  analysis  studies.  Miller 
(Ref.  1 )  and  Urick  (Ref.  2),  for  example,  proposed  theoretical  models  which  require  a  zero 
sound-speed  gradient.  These  models  cannot  account  for  surface-generated  noise  arriving  at 
angles  below  the  horizontal  or  for  bottom  noise  (biological,  seismic)  arriving  at  angles  above 
the  horizontal,  except  by  reflection  from  the  bottom  or  surface,  respectively.  This  results 
in  a  distorted  directional  noise  pattern  with  abnormally  low  levels  near  the  horizontal. 
Talham’s  model  (Ref.  3),  on  the  other  hand,  includes  realistic  sound-speed  profiles  but 
applies  only  to  bottom-mounted  hydrophones.  Bartberger  (Ref.  4)  considered  only  ship 
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l-igure  2  Typical  ambient  noise  vertical  directionality  plot. 


noise  from  a  uniform  distribution  of  surface  ships,  thus  restricting  the  usefulness  of  his 
model  to  the  20-  to  100-Hz  region.  Finally,  the  MOSS  Ambient  Noise  Model  developed 
jointly  by  the  Navy  and  Bell  Telephone  Laboratories  calculates  the  expected  vertical 
directionality  of  the  noise  field  for  seven  ocean  areas  for  the  mean  summer  and  winter 
conditions  for  acoustic  propagation,  wind  speed,  and  shipping.  A  more  general  area- 
independent  model  is  desired. 

In  addition  to  the  limitations  already  discussed,  none  of  the  above  models  considers 
horizontal  directionality,  noise  due  to  distant  sources  (SOFAR  Channel  noise),  or  noise  of  a 
biological  or  transient  nature,  nor  does  any  contain  sufficient  provision  for  considering  vari¬ 
ances  in  system  characteristics.  For  these  reasons  RAND1  was  developed  for  use  in  the  de¬ 
sign  of  ambient  noise  measurement  experiments  and  the  analysis  of  surveillance  systems. 


AMBIENT  NOISE  DIRECTIONALITY  MODEL 
MODEL  DESCRIPTION 

RANDI  utilizes  one  or  two  of  three  different  sources  for  the  propagation  loss  between 
the  noise  generators  and  the  sensor.  The  first  source  is  a  self-contained  linear  raytrace  routine, 
one  which  approximates  the  sound-speed  profile  by  a  series  of  straight-line  segments  and  cor¬ 
rects  for  earth  curvature.  See  Ref.  5  for  a  detailed  description. 

The  second  source  is  a  set  of  propagation  loss  versus  range  and  arrival  angle  arrays 
which  is  input  as  data.  If  the  propagation  loss  is  input,  RANDI  will  bypass  the  raytrace 
routine.  Hence,  ambient  noise  calculations  by  RANDI  can  be  based  on  propagation  loss 
values  from  ray  theory  (its  own  or  nearly  any  other  raytrace  model,  including  those  that 
account  for  variable  bottom  topography  and  horizontal  changes  in  the  sound-speed  profile), 
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normal-mode  theory,  experimental  measurements,  or  any  method  by  which  propagation 
loss  versus  range  and  arrival  angle  might  be  obtained.  This  second  source  of  propagation 
loss  can  be  extremely  useful  when  operating  RANDI  at  the  relatively  low  frequencies  of 
many  passive  surveillance  systems,  where  the  validity  of  ray  theory  becomes  questionable 
and  normal-mode  theory  more  attractive. 

The  third  form  of  propagation  loss  utilized  by  RANDI  is  for  sound  energy  traveling 
from  distant  sources  to  a  sensor  located  within  the  SOFAR  Channel.  RANDI  calculates  the 
propagation  loss  in  this  case  by  considering  the  effects  of  frequency-dependent  attenuation 
and  spreading  loss  increasing  with  fifteen  times  the  log  of  the  average  range  to  the  conti¬ 
nental  shelf  or  the  range  to  where  the  SOFAR  Channel  reaches  the  surface  (as  in  northern 
latitudes). 

In  addition  to  distant  noise,  RANDI  considers  five  other  sources  of  isotropic  and 
anisotropic  surface  and  volumetric  noise  (Fig.  3):  shipping,  sea  state  0,  biological,  rain,  and 
wind-wave  interaction. 

The  surface  noise  is  generated  by  an  infinite  number  of  point  sources  distributed 
along  a  horizontal  “noise  source  plane”  just  below  the  surface  (Fig.  4).  The  depth  of  this 
plane  is  set  at  20  ft,  since  the  main  ship  noise  radiators  (screws,  shaft,  and  hull)  are  near 
this  depth,  and  surface  wave  action  extends  well  below  the  surface.  The  shipping  noise 
generators  are  nonuniformly  distributed,  while  the  wind-wave  and  rain  noise  generators  are 
uniformly  distributed.  The  noise  resulting  from  any  of  these  three  sources  is  anisotropic. 

Sea  state  0  and  biological  noises  come  from  a  uniform  volumetric  distribution  of  noise  gen¬ 
erators  centered  around  the  sensor  and  result  in  isotropic  noise.  The  SOFAR  Channel  noise 
received  by  a  sensor  located  within  the  SOFAR  Channel  is  also  isotropic  for  all  channeled 
ray  angles.  The  total  noise  field,  then,  is  part  isotropic  and  part  anisotropic. 

A  target  capability  is  also  included  in  RANDI.  By  specifying  the  oceanographic, 
environmental  and  noise  conditions,  and  a  target  location,  depth,  and  frequency  spectrum 
or  line  component,  the  Surveillance  System  Analyst  receives,  by  way  of  plot  and  printed 
output,  the  levels  and  angles  of  target  multiple  arrivals  superimposed  on  the  ambient  noise 
arrivals  (see  Figs.  B-4  through  B-13).  Such  a  capability  can  make  RANDI  a  useful  tool  in 
the  design,  analysis,  and  optimization  of  surveillance  systems. 

MATHEMATICAL  DESCRIPTION 

The  anisotropic  shipping  noise  squared  pressure  spectrum  level  SNL  (for  a  1-Hz  band) 
in  the  model  is  obtained  from  surface  noise  generators  which  have  a  radiated  sound  pressure 
level  varying  with  frequency.  The  following  empirical  expression,  which  has  characteristics 
similar  to  the  spectra  given  in  Fig. .71  of  Ref.  6,  is  utilized: 

SNL  (re pPa)  —  A0  —  1 0  log  (TO-10  lQg f'+i  ]6  +  io+3.3  iogf-6.27) 

+  4  SHIPD  +  50  log  (SPEED/ 1 2)  +  20  log  (LENGTH/300) 
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where 

_  |  90  no  distant  noise  sources 
0  1  85  distant  noise  sources  specified 

f '  -  f  -  2  SHIPD  +  1  2 
f  =  number  of  hertz 

SHIPD  =  shipping  density  indicator  (0  —  7) 

SPEED  =  number  of  knots  of  noise-source  ships’  average  speed 
LENGTH  =  number  of  feet  of  noise-source  ships’  average  length 

The  squared  pressures  obtained  from  the  levels  given  by  this  function  are  distributed  over  the 
noise  source  plane  in  such  a  way  to  result. in  effective  squared  pressures  for  unit  area  at  unit 
distance  (the  mechanics  are  described  in  the  section  on  model  calibration)  with  magnitude 
varying  in  range  and  azimuth  in  such  a  manner  to  be  proportional  to  the  density  of  an  input 
shipping  density  distribution  ASHIP  (ship  density  vs  range  array)  or  a  shipping  density  in¬ 
dicator  SHIPD  (no  shipping  to  heavy  shipping  on  a  scale  from  zero  to  seven).  A  similar 
expression  was  used  for  the  squared  pressure  level  due  to  distant  sources  DSL,  since  ship¬ 
ping  is  the  major  contributor  at  low  frequencies  and  the  higher  frequencies  are  severely 
attenuated  at  SOFAR  Channel  propagation  ranges. 

The  squared  pressure  spectrum  levels  for  the  anisotropic  noise  WNL  due  to 
wind -wave  interaction  were  obtained  from  a  fit  to  the  data  of  Perrone  (Ref.  7)  and  are  a 
function  of  wind  speed  and  frequency  as  follows: 
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WNL  (re /uPa)  =  ~18  logf  +  96 


+  [9.66  (logf)2  -26  logf  +  27.3](0.065  WSPD07) 
where  WSPD  =  number  of  knots  of  wind  speed  in  the  vicinity  of  the  sensor. 

These  levels  are  converted  into  mean  squared  sound  pressure  in  a  1-Hz  band  and 
further  modified  by  a  multiplying  factor  to  achieve  directionality  of  the  source  at  high  fre¬ 
quencies  and  nondirectionality  at  the  lower  frequencies  as  indicated  by  Becken  (Ref.  8). 
The  wind  noise  squared  pressure  modifier  used  in  the  model  is  the  following: 

Squared  pressure  multiplier  =  1  —  (1  —  cosn  0)  (0.002  f  —  0.02) 


where  0  =  angle  (in  degrees)  from  the  vertical  that  the  ray  makes  at  the  source,  and 


90-0 
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0  <  0  <  70 
70  <  0  <  90 


The  following  function,  obtained  from  the  data  of  Franz  (Ref.  9),  gives,  for  rain, 
the  expected  noise  squared  pressure  level  RPL  in  a  1-Hz  band  and  is  distributed  uniformly 
along  the  horizontal  noise  plane  in  a  manner  similar  to  the  wind-wave  noise. 

RPL  (re  MPa)  -  5.5  logf  +  50.5  +  14.5  log  (RAIN) 
where  RAIN  is  the  number  of  inches  of  rainfall  per  hour  in  the  vicinity  of  the  sensor. 

Isotropic  sea  state  0  noise  squared  pressure  spectrum  level  SS0  is  that  noise  which  is 
measured  under  the  ideal  conditions  of  no  wind,  calm  surface,  no  biological  activity,  and 
negligible  shipping.  It  varies  with  frequency  and  is  independent  of  depth  and  geographic 
location.  The  equation  for  SS0  (for  a  1-Hz  band)  was  obtained  by  a  quadratic  fit  to  data 
of  Wenz  (Ref.  10): 

SS0  (re  pPa)  -  4.22  (log  f)2  -  33.4  logf  +  89.1 

A  volumetric  or  isotropic  distribution  was  chosen  for  biological  noise  in  the  absence 
of  data  indicating  otherwise.  This  noise  varies  with  time  of  day  and  relative  amount  of  ac¬ 
tivity  expected  at  a  particular  site.  An  activity  indicator  on  a  scale  from  zero  to  ten  is  an 
input  to  the  model.  The  equation  for  the  biological  noise  squared  pressure  spectrum  level 
BPL  is  of  the  following  form: 

BPL  (re  MPa)  -  0.00175  (100  -  f)  sin  [0.00262  (hr- 300)]  ACTIVITY 
-  16.96  (logf)2  +  50.1  logf  +  45.1 
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where  the  0<  ACTIVITY  <  10  and  hr  is  the  local  time  of  day  (military  designation 
0000-2400). 

The  sound  pressure  spectrum  level  TSL  for  a  target  source  is  accepted  by  the  model 
by  reading  in  the  appropriate  coefficients  A(,  i  =  1 , 3,  of  a  third-degree  polynomial 

TSL  =  A3  (log  f)3  +  A,  (log  f)2  +  Aj  log  f  +  A0 

In  the  event  that  the  tine  component  is  the  dominant  feature  in  the  band  of  interest,  A0  is 
the  level  of  the  line.  Also  specified  by  the  user  are  the  initial  target  range  and  depth,  the 
numbers  of  target  ranges  and  depths  to  be  calculated,  and  the  increments  in  range  and  depth 
to  be  used. 

COMPUTATIONAL  PROCEDURE 

All  of  the  noise  and  target  squared  pressure  spectrum  levels  are  integrated  over  a 
user-specified  bandwidth  by  means  of  an  input  frequency,  response  function.  If  none  is 
specified,  a  1-Hz  bandwidth  and  a  constant  bandpass  frequency  response  function  are 
assumed. 

The  squared  pressure  received  for  a  differential  vertical  angle  is  obtained  by  first 
calculating  the  area  defined  by  the  intersections  of  the  corresponding  ray  bundle  with  the 
noise  source  plane.  These  areas  are  multiplied  by  the  local  effective  squared  noise  pressures 
for  unit  area  at  a  distance  of  1  yard.  The  resulting  squared  pressures  are  then  reduced  to 
account  for  propagation  loss  and  summed.  To  this  value  is  added  the  contribution  of  the 
isotropic  noise  sources.  The  result  is  further  reduced  to  account  for  the  vertical  response  of 
an  individual  hydrophone  or  an  array  of  hydrophones,  If  no  response  function  is  specified, 
it  is  assumed  omnidirectional.  The  final  squared  noise  pressure  arriving  at  that  angle  is 
stored  for  output  and  the  process  is  repeated  for  an  adjacent  ray  bundle  at  a  new  vertical 
angle. 


In  the  case  where  horizontal  directionality  is  desired,  the  ocean  is  divided  into  n 
regions  by  passing  vertical  planes  through  the  sensor  location  at  360/n  deg.  The  ocean  is 
then  effectively  divided  into  n  regions,  similar  to  a  huge  sliced  pie  of  infinite  radius  whose 
thickness  is  equal  to  the  ocean  depth  at  the  receiver  location.  The  pertinent  environmental 
and  noise  parameters  are  specified  separately  for  each  “pie  slice”  region.  Hence,  the  calcu¬ 
lations  performed  for  one  region  are  independent  of  the  calculations  performed  in  adjacent 
regions.  The  total  squared  noise  pressure,  that  which  would  be  measured  by  an  omnidi¬ 
rectional  hydrophone,  is  obtained  by  summing  the  squared  noise  pressures  for  the  n  inde¬ 
pendent  “pie  slice”  regions.  An  example  is  given  in  Appendix  B,  Figs.  B-6  through  B-l  1 . 

An  explanation  of  terminology  used  in  the  ambient  noise  directionality  illustrations 
is  necessary  before  results  can  be  interpreted.  Noise  level  refers  to  the  per  steradian  mean 
squared  pressure  spectrum  level.  Vertical  received  angle  is  the  angle  at  the  sensor  the  inci¬ 
dent  ray  makes  with  the  horizontal.  The  negative  rays  are  downcoming  rays  at  the  sensor, 
with  the  ray  having  a  vertical  received  angle  of  -90  deg  being  the  ray  which  arrives  from 
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the  surface  directly  above  the  sensor.  The  rays  with  positive  angles  arrive  at  the  sensor 
from  angles  below  the  horizontal. 

Unlike  noise,  targets  are  treated  as  point  sources  at  specific  ranges  and  depths.  How¬ 
ever,  the  received  signal  squared  pressures  are  similarly  influenced  by  the  bandwidth,  fre¬ 
quency  response  function,  and  the  vertical  beam  response  function. 

To  aid  the  user,  many  of  the  model  inputs  have  been  initialized  in  data  statements. 
This  eliminates  the  repeated  input  of  variables  which  are  common  to  many  situations  and 
yet  allows  the  initial  values  to  be  suppressed  in  the  event  that  different  values  are  desired. 
Those  parameters  which  have  been  initialized  are  given  in  the  model  input  section  of  the 
program  listing,  Appendix  C. 

MODEL  CALIBRATION  PROCEDURE 

The  functional  relationships  used  for  noise  express  the  pressure  spectrum  levels  one 
would  expect  to  measure  and  are  not  source  levels.  To  get  the  effective  source  levels  at  a 
distance  of  1  yard  (for  unit  surface  area)  for  distributions  of  noise  generators  requires  re¬ 
moving  from  the  original  levels  the  effects  of  frequency-dependent  attenuation  in  the  me¬ 
dium  and  then  normalizing  to  a  unit  of  surface  area.  This  is  accomplished  by  distributing 
the  noise  generators  along  the  horizontal  noise  plane  and  propagating  the  noise  at  different 
frequencies.  A  calibration  function  can  then  be  obtained  which  is  added  to  the  original 
levels  to  yield  corrected  levels  and  eliminate  bias.  This  effectively  removes  the  effects  of 
frequency  attenuation  in  the  original  levels  and  converts  them  to  source  levels  for  unit  sur¬ 
face  area.  Finally,  the  output  of  the  model  is  adjusted  to  coincide  with  measured  data  for 
one  set  of  conditions  (frequency,  depth,  wind  speed,  etc.)  at  one  location. 

This  calibration  is  dependent  on  the  manner  in  which  the  noise  sources  arc  distributed 
and  the  method  by  which  the  noise  generation  area  is  calculated  and  is  independent  of  the 
medium.  Hence,  once  this  calibration  is  performed  for  one  set  of  conditions  at  one  location, 
it  need  not  be  done  for  other  conditions  or  locations.  The  initial  calibration  was  performed 
using  Marine  Physical  Laboratory  (MPL)  Pacific  ambient  noise  data  (Ref.  1  1)  and  shipping 
information  provided  by  Western  Seas  Frontier. 

MODEL  VALIDATION 

The  final  step,  validation,  has  not  yet  been  completed.  This  requires  comparing  the 
model  output  with  measured  data  for  many  different  locations,  seasons,  depths,  frequencies, 
shipping  distributions,  wind  speeds,  etc.  Tentative  plans  include  comparing  with  the 
IOMF.DEX  data,  which  were  taken  in  the  Mediterranean  during  November  1971 .  Com¬ 
parison  with  other  data  sets  for  different  conditions  will  be  done  as  data  become  available. 
Although  a  large  quantity  of  data  is  now  readily  available,  the  lack  of  shipping  information 
for  the  time  interval  during  which  the  measurements  were  taken  disqualifies  it.  Nearly  all 
historical  ambient  noise  data  are  inappropriate  for  model  validation. 
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DISCUSSION  OF  RESULTS 


GENERAL  RESULTS 

Examples  of  output  from  the  model  for  three  different  geographic  locations  (Fig.  5) 
have  been  included  to  illustrate  both  the  capability  of  the  model  and  variability  in  the  di¬ 
rectional  character  of  the  noise  field  resulting  from  differences  in  sensor  depth  and  acoustic 
propagation  (Figs.  6,  7,  and  8).  Profiles  which  generally  characterize  summer  conditions  in 
the  Pacific.  Atlantic,  and  Mediterranean  were  chosen. 

It  is  interesting  to  note  the  existence  of  the  trough  (low  level  of  noise)  in  the  noise 
field  near  the  horizontal.  This  condition  results  from  the  horizontal  rays  at  the  hydrophone 
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being  unable  to  reach  the  surface.,  where  the  wind  and  shipping  noise  is  generated,  In  two 
cases,  as  the  depth  increases,  the  trough  becomes  narrower  and  eventually  disappears,  indi¬ 
cating  that  only  near  the  bottom  will  the  horizontal  rays  “see”  the  surface.  In  the  other 
case,  however,  the  trough  gets  wider  as  the  depth  increases.  This  results  from  the  sound 
speed  at  the  bottom  being  much  less  than  that  at  the  surface,  hence,  horizontal  bottom  rays 
can  never  reach  the  surface.  Energy  within  the  trough  arrives  from  distant  sources.  The 
trough  levels  are  determined  by  the  amount  of  energy  which  has  become  trapped  within  the 
SOFAR  Channel.  There  can  be  a  trough,  or  it  can  be  filled  in,  depending  on  whether  the 
energy  per  arrival  angle  due  to  near  sources  is  greater  or  less  than  that  due  to  distant  sources. 
The  total  amount  of  sound  energy  received  from  the  distant  sources  depends  on  the  total 
number  of  ships  over  the  continental  slope  and  shelf,  range  to  the  shelf,  and  the  total  channel 
“look”  angle  at  the  sensor  depth. 

The  lower  level  of  sound  at  down  going  angles,  greater  than  approximately  1  5  deg, 
compared  to  the  up  going  angles,  less  than  —15  deg,  results  from  bottom  loss.  The  upcoming 
rays,  having  bounced  off  the  bottom  before  reaching  the  receiver,  suffer  a  loss  in  energy  not 
experienced  by  rays  coming  from  the  surface.  This  accounts  for  the  asymmetric  shape  of 
the  noise  plots. 

MODEL  OUTPUT-MEASURED  DATA  COMPARISON 

Ambient  noise  data  obtained  by  MPL  (Ref.  1 1 )  during  their  April  and  May  1 97 1 
FLIP  ambient  noise  measurements  were  utilized  for  the  initial  model  calibration  and  to  per¬ 
form  a  level-variable  dependence  check. 

The  model  was  initially  calibrated  at  50  Hz  and  100  m  receiver  depth  by  reading  in 
the  measured  wind  speed  and  recorded  shipping  data  and  adjusting  the  model  output  to 
correspond  to  the  measured  noise  level  for  a  single  day  of  the  15-day  series  of  measure¬ 
ments.  The  output  of  the  model  was  then  compared  with  the  measured  data  for  other  fre¬ 
quencies,  depths,  wind  speeds,  and  shipping  distributions. 

Figures  9a  and  b  are  examples  of  bow  RANDI  compares  with  the  MPL  data  as  a 
function  of  time.  Only  every  other  point  of  the  original  MPL  data  (one  each  day)  has  been 
retained  to  compare  with  the  model,  since  the  shipping  distribution  near  the  FLIP  site  was 
reported  only  once  each  day.  Figure  9a  is  a  comparison  of  the  MPL  data  and  the  RANDI 
output  for  200  Hz  at  560  m  receiver  depth.  This  frequency  is  of  particular  interest  in 
checking  the  output  of  the  model,  since,  at  this  frequency,  both  shipping  and  wind-wave 
interaction  make'  significant  contributions  to  the  total  noise  level.  This  figure  indicates  that 
the  model  is,  in  general,  within  2  or  3  dB  of  the  measured  data.  The  anomalous  results  can 
be  explained  by  the  presence  or  lack  of  nearby  (within  the  same  1-deg  square  as  FLIP)  ship¬ 
ping  which  was  reported  (erroneously  perhaps)  but  was  not  seen  by  FLIP  personnel  or  could 
have  been  close  to  FLIP’S  position  but  not  reported.  The  agreement  at  400  Hz  (Fig.  9b)  is 
similarly  encouraging  for  such  a  rough  initial  calibration. 

Figure  9c  compares  the  noise  levels  calculated  for  MPL  by  CAPT  Paul  Wolff 
(Ref.  11)  using  two  different  noise  models,  the  noise  as  calculated  by  RANDI  and  the 
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Figure  9.  RANDI  model-measured  data  comparison  for  several  frequencies  and  depths  (Ref.  1 1). 


MPL  data  for  100  Hz  at  100  m  depth,  Again  the  agreement  between  RAND1  and  the 
MPL  data  is  good. 

RANDI  also  indicated  that  the  noise  field  should  have  a  depth  dependence.  The 
total  reduction  in  noise  level  as  the  sensor  depth  increased  from  100  to  3930  m  was  calcu¬ 
lated  to  be  6  dB  at  50  Hz  and  4  dB  at  400  Hz.'  The  data,  however,  showed  a  reduction  of 
about  10  and  4  dB,  respectively,  for  the  same  conditions.  The  discrepancy  at  50  Hz  is 
largely  attributable  to  not  specifying  distant  (SOFAR  Channel)  noise  sources  in  RANDI 
when  the  comparisons  were  made.  This  could  easily  add  the  required  4  dB,  since  the 
total  SOFAR  Channel  “look  angle”  at  the  axis  is  approximately  22  deg  and  the  nearby 
shipping  was  relatively  light.  A  similar  depth  dependence  was  observed  by  Lomask  and 
Frassetto  in  the  Mediterranean,  where  they  utilized  the  bathyscaph  Trieste  (Ref.  1  2),  Weigle 
and  Watt  (Ref.  13),  however,  generally  observed  less  than  3  dB  depth  dependence  in  the 
Mediterranean  for  the  same  frequency  range.  The  environmental  and  shipping  data  from 
Ref.  13  were  used  as  input  for  RANDI  to  enable  comparison  with  the  measured  data.  The 
results  are  presented  in  Table  1 .  Since  the  data  are  classified,  only  the  differences  in  the 
measured  and  calculated  (by  RANDI)  levels  are  given  as  a  function  of  frequency  and  depth. 
It  is  interesting  to  note  that  although  RANDI  was  calibrated  with  Pacific  data,  where  the 
ambient  noise  levels  are  about  20  dB  less  than  in  the  Mediterranean,  RANDI  calculated  the 
Mediterranean  noise  within  about  1  or  2  dB. 

RANDI  was  also  used  to  obtain  the  horizontal  directionality  of  the  ambient  noise 
for  one  shipping  distribution  and  wind  speed  during  the  FLIP  ambient  noise  measurement 
program  The  results  are  given  in  Figs.  B-6  through  B-l  1 . 


FUTURE  IMPROVEMENTS 

Improvements  in  the  current  ambient  noise  directionality  model  will  be  made  as 
time  and  state-of-the-art  permit. 

Transient  noise,  including  noise  from  seismic  disturbances  and  explosions,  is  of  con¬ 
siderable  interest  and  will  be  included  in  RANDI  within  the  immediate  future. 


SUMMARY 

RANDI  is  a  working  FORTRAN  model  which  calculates  and  plots  the  directional 
character  of  low-frequency  ambient  noise  in  the  ocean  environment.  This  model  is  con¬ 
ceptually  simple  and  will  be  improved  as  time  and  state-of-the-art  permit.  It  is  intended  to 
be  a  tool  to  aid  in  the  design  of  noise  measurement  experiments  and  in  surveillance  systems 
design  and  evaluation  and  not  for  ambient  noise  synoptic  forecasting. 
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Table  1.  Comparison  oF  RAND1  Output  and  Measured 
Noise  Pressure  Spectrum  Level  Data  Prom  the 
Mediterranean  (Ref.  13). 


Frequency,  Hz 

Depth,  ft 

Difference  in  Level,  dB 

20 

925 

-1.0 

1728 

-0.6 

3785 

5795 

-1.9 

6730 

-1.8 

SO 

925 

1.9 

1728 

2.6 

3785 

0.4 

5795 

2.7 

6730 

4.5 

100 

925 

0.5 

1728 

0.0 

3785 

0.3 

5795 

1,1 

6730 

2.4 

200  ■ 

925  ■ 

-1.2 

1728 

-1.6 

3785 

-2.8 

5795 

-0.5 

6730 

-0.1 

400 

925 

02 

1728 

0.3 

3785 

-1.6 

5795 

-1.6 

6730 

-0.4 

500 

925 

0.5 

1728 

-0.7  ■ 

3785 

-0  5 

5795 

-16 

6730 

-11 
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Appendix  A 


RANDI  SUBROUTINE  DESCRIPTIONS  AND  FLOW  DIAGRAM 


Figure  A-l  illustrates  how  the  main  or  executive  routine  controls  the  flow  of  logic 
through  the  21  subroutines  which  comprise  RANDI.  This  figure  can  also  be  used  as  a  guide 
in  '‘streamlining”  RANDI  by  removing  those  subroutines  which  will  not  be  required  by  the 
user.  In  such  a  case,  “dummy”  or  “empty”  subroutines  are  substituted  for  the  originals. 

For  example,  RANDI  could  easily  be  loaded  into  a  relatively  small  computer  when  the  prop¬ 
agation  loss  versus  range  is  to  be  read  in  and  no  sound-speed  profile  plot  is  required.  This 
is  accomplished  by  substituting  “empty”  subroutines  for  TARGET,  RAPATH,  SUMINT, 
RAYTRC,  VERTEX,  HSURF,  SCLNCE  and  SSPLOT.  If  only  a  printout  is  desired,  PLT 
could  also  be  deleted. 


The  following  is  a  brief  functional  description  of  each  subroutine: 


ALFA  contains  the  equation  for  the  absorption  coefficient  as  a  function  of  frequency 
and  water  temperature. 

AUXPR  calculates  sea  surface  area  intersected  by  a  ray  bundle  and  sums  squared  sound 
pressure  when  propagation  loss  is  read  in. 


BIO 

BWIDTH 

ERTHC 

FUNS 

FUNU 

HLOSS 

HSURF 

PLT 

PRINTS 


contains  equation  for  biological  noise  as  a  function  of  frequency,  time  and 
biological  activity  indicator, 

integrates  noise  or  target  spectrum  using  an  input  bandwidth  response  function. 

corrects  the  sound  speed  profile  to  include  the  effects  of  earth  curvature. 

obtains  the  squared  noise  pressure  from  the  shipping  histogram  array  for  a 
given  range. 

linearly  interpolates  between  two  points. 

computes  absorption,  spreading  and  refraction  losses. 

computes  surface  reflection  losses. 

plots  noise  and  target  levels 

prints  input  and  initial  variables. 
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READIN 


Figure  A-l .  RANDI  subroutine  How 


RAPATH 

RAYTRC 

REA  DIN 

SCLNCE 

SHIPIN 

SSPLOT 

SUMINT 

SURFI 

TARGET 

VERTEX 


generates  ray  paths  and  corresponding  propagation  losses, 
determines  path  taken  by  each  ray. 
reads  input  parameters. 

produces  convenient  set  of  scale  numbers  for  sound-speed  profile  plot, 
constructs  the  ray  for  the  shipping  noise  squared  pressure  versus  range  histogram, 
plots  sound-speed  profile. 

calculates  sea  surface  area  intersected  by  a  ray  bundle  and  sums  squared  sound 
pressure  when  propagation  loss  not  read  in. 

contains  the  squared  pressure  levels  for  thermal,  shipping,  wind-wave  and  rain 
noise  generators. 

traces  rays  from  target  to  sensor  and  computes  received  signal  levels  and  arrival 
angles. 

computes  greatest  and  shallowest  depths  reached  by  ray 
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Appendix  B 

RAND1  INPUT-OUTPUT 


RANDI  INPUT  DATA  AND  FORMAT 


Input  data  cards  for  RANDI  include  parameter  cards,  array  cards  and  run-control 

cards. 

A  RANDI  input  parameter  card  consists  of  an  arbitrary  number  of  parameter  fields 
in  free  format.  Each  parameter  field  is  separated  by  a  comma,  oblique,  or  blank  space  and 
contains  the  parameter  name  followed  by  an  equal  sign  and  then  the  value  the  parameter  is 
to  be  assigned  for  that  and  subsequent  consecutive  runs  or  until  a  new  value  is  assigned.  The 
array  cards  contain  the  name  of  the  array  followed  by  an  equal  sign,  which  is  immediately 
followed  by  the  first  value  in  the  array.  Subsequent  values  are  separated  by  a  blank  space. 
Values  which  are  continued  on  additional  cards  are  preceded  by  an  asterisk  and  a  blank 
space  at  the  beginning  of  the  string  of  values  contained  on  the  continuation  card. 

Data  cards  for  consecutive  runs  are  separated  by  a  PAUSE  run  control  card,  which 
signals  the  end  of  data  for  one  run  and  that  further  data  follows  for  another  run.  The  last 
data  set  is  followed  by  an  END-DATA  card  signaling  the  end  of  data  for  the  last  run. 

The  PRLOSS  card  precedes  each  array  of  propagation  loss  versus  range  and  follows 
the  parameter  and  array  cards  for  a  given  run  but  precedes  the  PAUSE  card.  The  plot  data 
cards  follow  the  PAUSE  card  at  the  end  of  that  data  set  for  a  given  run,  since  they  are  read 
directly  by  the  plot  routines  and  not  the  READIN  subroutine. 


Although  it  is  possible  to  read  in  a  total  of  thirty-six  parameters  and  arrays,  RANDI 
can  be  run,  for  the  large  majority  of  cases,  by  merely  reading  in  a  few  variables.  For  ex¬ 
ample,  if  only  the  sound-speed  profile,  frequency,  and  sensor  depth  are  read  in,  RANDI 
will  use  preprogrammed  values  for  the  other  variables  and  arrays  in  calculating  the  ambient 
noise.  The  preprogrammed  values  are  listed  in  the  model  input  section  of  the  model  listing 
of  Appendix  C.  The  resulting  noise  calculation,  then,  will  be  for  the  given  depth,  a  wind 
speed  of  5  knots,  moderately  heavy  shipping,  an  omnidirectional  vertical  beam  response 
function,  and  a  1-Hz  bandwidth  centered  at  the  input  frequency  integrated  over  a  bandpass 
frequency  response  function.  A  360-deg  horizontal  sector  width,  moderate  bottom  loss, 
and  a  latitude  of  40-deg  will  also  be  used. 


The  following  is  a  list  and  brief  description  of  the  RANDI  input  variables  and  output 


INPUTS 


products. 

System 

Name 

Units 

ADANF 

ARRAY 

ZTCi 

FEET 

PHID 

DEG 

Frequency  (kHz)  and  depth  (ft)  array 
Noise  source  depth 
Declination-elevation  angle 
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Name  Units 

DELPH  DEG  Vertical  half-power  beamwidth 

DGREH  DEG  Width  of  horizontal  sector 

BERNG  DEG  Bearing  of  horizontal  sector 

BNCS  Number  of  ray  bounces 

ABW  ARRAY  Frequency  response  (dB  down)  across  bandwidth 

ARESP  ARRAY  Vertical  beam  response  (dB  down)  pattern 


ENVIRONMENTAL 

ALAT  DEG  Latitude  for  earth  curvature  effects 

AVELP  ARRAY  Sound-speed  (ft/sec)  vs  depth  (ft)  profile 

AHB  ARRAY  Bottom  reflection  loss  (dB)  vs  grazing  angle  (deg) 

APROP  ARRAY  Optional  prop  loss  (dB)  vs  range  (kyd)  with  elevation  angle  (deg) 

and  ray  bundle  width  (deg) 


NOISE  SOURCES 


WSPD 

KT 

Wind  speed  in  knots 

ACTIV 

Biological  activity 

HOUR 

HR 

Local  time  of  day  in  hours  (military) 

RAIN 

IN/HR 

Rainfall 

SHIPD 

Ship  density  scale 

ASHIP 

ARRAY 

Ship  density  (No. /1 0000  sq  mi)  vs  range  (nmi) 

SHLFR 

NM 

Average  distance  to  continental  shelf  and  surfaced  sound  channel 

SHLFS 

Number  of  ships  within  horizontal  sector  over  continental  shelf 
and  in  surfaced  sound  channel 

TARGET 


ZTG1 

FT 

Target  depth 

ZTGNO 

Number  of  target  depths 

ZTING 

FT 

Target  depth  increment 

RGT 

KYD 

Range  to  target 

RGTNO 

Number  of  target  ranges 

RGINC 

KYD 

Target  range  increment 

TARG 

If  less  than  0  TARGO  is  a  target  line  component  in  dB/micropascal, 
otherwise  TARGO  is  const  in  target  spectrum  level  equation 

TARGO 

Level  of  line  component  in  dB  or  const  in  spectrum  equation 

TARG1 

Coefficient  of  log  f 

TARG2 

Coefficient  of  (log  f)2 

TARG3 

Coefficient  of  (log  f)3 

24 


OUTPUT  AND  PLOT  CONTROL  CARDS 


System 

Name 

OUTPT 

Output  data  control  parameter 

SNPLT 

Noise  plot  flag 

HFLAG 

Noise  plot  data  card  flag 

SSPLT 

Sound-speed  profile  plot  flag 

Three  additional  cards  required  after  pause  for  noise  plot 
Title 
Location 
Date  . 

Two  additional  cards  required  for  sound-speed  plot 
Location 
Date 


INPUT  DATA  CONTROL  CARDS 


HEADER 

PAUSE 

PRLOSS 

END-DATA 


Used  before  a  message  statement 

Used  before  each  consecutive  run  except  first 

Used  before  each  prop  loss  array 

Used  at  end  of  data  or  before  plot  cards  that  are  at  end  of  data 


OUTPUTS 


Plots 

1.  Noise  level  versus  vertical  arrival  angle  with  or  without  target 

2.  Sound-speed  profile 

Table 

Noise  level  versus  vertical  arrival  angle 
Target  signal  level  versus  vertical  arrival  angle 
Total  sector  noise  level 


EXAMPLE  COMPUTER  RUN 

Input 

An  example  of  a  typical  RANDI  run  follows  to  illustrate  the  input  format  and  out¬ 
put  products.  In  this  example  both  a  vertical  directionality  plot  and  a  sound-speed  profile 
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plot  are  requested.  The  input  variables  are  sensor  depth,  frequency,  latitude,  wind  speed, 
request  for  both  plots,  bottom  loss  versus  grazing  angle  array  and  sound  speed  versus  depth 
array.  Three  plot  data  cards  are  also  included,  giving  the  title,  location,  and  date.  All  other 
parameters  are  automatically  set  equal  to  preprogrammed  initial  values. 

Output 

Figures  B-2  through  B-5  illustrate  the  output  of  RANDI  for  a  given  set  of  input 
data  cards. 


EXAMPLE  OF  HORIZONTAL  DIRECTIONALITY 

Shipping  data  for  May  1 .  1  971 ,  within  a  few  hundred  miles  of  the  FLIP  ambient 
noise  measurement  site,  were  utilized  to  illustrate,  the  horizontal  directionality  capability 
of  RANDI.  To  do  this  the  noise  field  was  divided  into  36  horizontal  pie-shaped  sectors 
centered  at  the  FLIP  position.  The  vertical  arrival  structure  of  the  noise  field  was  obtained 
by  running  RANDI  with  different  shipping  distributions  each  time  to  get  the  average  noise 
level  every  10  deg.  These  levels  were  then  reduced  by  10  dB  to  get  the  per  degree  levels 
and  plotted  by  a  separate  DISSPLA  plot  routine.  Four  of  the  1 0-deg-horizontal-sector 
vertical  directionality  plots  (Figs.  B-8-1 1 )  have  been  included  to  demonstrate  the  variability 
in  the  vertical  arrival  structure. 
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Figure  B-l-  Example  run  input  data  cards. 
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Figure  B-2.  Typical  RANDI  initial  parameter  printout. 


LOCATION  06l00w  2400N 
DATE  03  23  72 
WIND  SPEED  (KT)  5.0 
BOTTOM  DEPTH  (FT)  19685.0 


FREQ  (HZ)  50.0 
REC  DEPTH  (FT)  975.0 
SECTOR  LEVEL  (DB)  80.3 
HOR  SEC  (DEG)  360  AT  270 
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Figure  B-3.  Typical  RANDI  noise  field  vertical  directionality  printout. 


NOISE  LEVEL  CDB  RE  mPH/STR  HZ) 


EXAMPLE  NOISE  PLOT 


LOCUTION  061 OOM  2«QN  FREQ  (HZ)  50.0 
06IE  03  23  72  RED  DEPTH  HU  975.0 
HIM)  SPEED  OCT)  5.0  SECTOR  LEVEL  (OB)  80.3 


HORIZONTRL  DIRECTIONALITY 

KPL  PflClFIC  DflTn  -  5/1/71 
100  HZ  FIT  300  n,  25  KT  RIND 
80.62  DB  OMNI  LEVEL 


180 


Figure  B-6.  Pei  degree  horizontal  directionality  of  ambient  noise  in  1  0-deg  sectors, 
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latitude  N,  degrees 


longitude  W,  degrees 


32 

30 

28 

26 

24 


127  125  123  121  119 


r...  . — -T- - T" 

315"  15"  45' 


figure  B-7.  Pacific  shipping  data  for  1  May  1971. 
Number  of  ships  over  1000  gross  tons  per  J-deg 
square. 
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Figure  H-t>.  Ambient  noise  vertical  directionality  for  the  10-deg  horizontal  sector  at  15  deg. 


NOISE  LEVEL  CDB  RE  MPA  /SIR  HZ) 


MPL  DRTR  COMPARISON 


LOCATION  2800N  12300U 
DATE  05  01  71 
HIND  SPEED  (NT)  25.0 
BCTTOn  DEPTH  <PT)  13615.0 
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SECTDK  LEVEL  1DB)  66.2 
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Figure  B-9.  Ambient  noise  vertical  directionality  lor  the  10-deg  horizontal  sector  at  45  deg. 


NOISE  LEVEL  CDB  RE  (jPfl  /STR  HZ] 


MPL  DATA  COMPARISON 


ON 


LOCATION  2800N  12300H  PREO  (HZ)  100.0 
DATE  050171  REC  DEPTH  (FT)  300.0 
HIND  SPEED  1KT1  25.0  SECTOR  LEVEL  (DB)  46.5 
BOTTOM  DEPTH  (FT)  13615.0  HOR  SEC  (0E6)  10  AT  12S 
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Figure  B-10.  Ambient  noise  vertical  directionality  for  the  10-deg  horizontal  sectorat  125  deg. 
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^  Figure  B-ll.  Ambient  noise  vertical  directionality  for  the  1 0-deg  horizontal  sector  at  315  deg. 
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TARGET 


LOCHTION  06100H  2400N 
DHTE  03  23  72 
HIM)  SPEED  (KT1  5.0 

Borron  depth  tra  19005.0 


Figure  B-l 2.  Example 


PLOT  EXAMPLE 


FREQ  (HZ)  50.0 
REC  OEPW  (FT)  875.0 
SECTOR  LEVEL  (DBJ  80.3 
HOR  SEC  (DEB)  360  RT  270 


TARGET  DATA 

-  400  FT  TURBO  HI  50.0  KTD 

58.1  06  REC  FROM  150.0  D0  SOU 
52.3  00  HRX  ONE  DECREE  flVERHOE 


^ED  ANGLE  (DEG) 

et— ambient  noise  plot. 
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Figure  B-13.  Example  target  data  output. 
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Appendix  C 


RESEARCH  AMBIENT  NOISE  DIRECTIONALITY 
(RANDI)  MODEL  LISTING 


oonnooooo 


c  *********************************************************** 

C  *  * 

c  *  RESEARCH  AMBIFNT  NOISE  DIRECTIONALITY  (RANDI )  MODEL  * 


C  ****************  **********************************  ********* 


c 

c 

c 

c 

*****  model  INPUTS  ***** 

c 

c 

c 

name 

UNITS 

initial 

C** 

c 

SYSTEM  INPUTS 

** 

c 

c 

adanf 

ARRAY 

N  FREQUENCIES  (KHZ)  AND  K  DEPTHS  (FT) 

c 

N+K  .LE.  30 

c 

ZTG 

feet 

NOISE  SOURCE  DEPTH  - 

20 

c 

PHID 

deg 

O/E  angle 

0 

c 

delph 

deg 

VERTICAL  HALF-POWER  BEAMW1DTH 

170 

c 

dgreh  ‘ 

DEg 

WIDTH  OF  HORIZONTAL  SECTOR 

360 

c 

BERNg 

DEg 

bearing  of  horizontal  sector 

270 

c 

BNCS 

BOUNCES  .LE .  60 

10 

c 

ABW 

array 

FREQUENCY  RESPONSE (DB  down)  across 

c 

FREQ  IN  BAND)  1  HZ  BW  INITIAL 

c 

aresp 

array 

VERTICAL  BEAM  RESPONSE  PATTERN 

c 

c 

OMNI  ASSUMED  IF  NOT  INPUT 

c 

c 

c** 

c 

environmental 

INPUTS  ** 

c 

c 

ALAT 

DEg 

latitude  for  earth  curvature  effects 

40 

c 

AVELP 

array 

VELOClTY(FT/SEC)  Vs  DEPTH(FEET)  profile 

c 

ANB 

array 

oOT-REFL  LOSS(DB)  VS  GRAZING  ANGLE (DEG) 

c 

PROV  5  ASSUMED  IF  NOT  INPUTTED 

c 

APROP 

array 

OPTIONAL  PROP  LOSS  (DB)  VS  RANGE 

c 

( K YD )  ARRAY  WITH  ELEVATION  ANGLE 

c 

(DEG)  AND  RAY  BUNDLE  WIDTH  (DEG) 

c 

APROP(l)  =  depression  angle  (DEG) 

c 

APROP(2)  =  RAT  BNDL  WDT  (DEG) 

c 

APROP ( 2N+1 )  =  RANGE  (KYDS) 

c 

APROP ( 2N+2 )  =  PROP  LOSS  (DB) 

C 


C**  NOISE  SOURCE  INPUTS  ** 


WSPD 

KT 

WIND  SPEED  IN  KNOTS 

5 

ACTIV 

BIOLOGICAL  ACTIVITY 

SCALE  0  TO  10 

0 

HOUR 

HR 

TIME  OF  DAY  IN  HOURS 

(MILITARY) 

1500 

RAIN 

IN/HR 

RAIN  FALL 

0 

SHI  PD 

SHIP  DENSITY  SCALE 

0  TO  7 

5.2 

.EG.  0  -  NO  SHIPS'  =  7  -  DENSE 

SHIPPING  CALIB  -  MPL  DATA  4/24/71 
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n  n  n  n  n 


ASHIP 


C 

C 

C 

c 


ARRAY  SHIPPING  HISTOGRAM  ARRAY  _  EACH  BAR  IS 

DEFINED  BY  2  To  4  NUMBERS  IN  THE 
FOLLOWING  ORDER  -  RANGE  IN  Nm> 

NO.  SHIPS  PAST  THAT  RANGE  AND 


C 

C 

C 

C 

c 

c 

c 

c 

c 

c 


before  next  range 

NEG  AV  SHIPS  SPEED  <KT>  t  OPTIONAL)  -12 
NEG  AV  SHIPS  LENGTH  (FT)  (OPT)  -300 

LAST  DATA  POINT  IS  LAST  RANGE 

SHLFR  NM  AVERAGE  DISTANCE  IN.  NAUTICAL  MILES  900 

to  the  continental  shelf  and 
surfaced  sound. CHANNEL  within  the 

SECTOR  of  HQRIZ  ANGLE  (DGREH) 

SHLFS  NO.  SHIPS  WITHIN  HORIZ  SECTOR  OVER  0 

CONT  SHELF  AND  IN  SURFACED  SOUND 
CHANNEL  (ATLANTIC  VLAM  DATA  CAL) 


C**  TARGET/THREAT  INPUTS  **  (WHEN  PROP  LOSS  NOT  READ  IN) 

C 


C 

ZTG1 

FT 

target  Depth  ,lt.  o  no  target 

c 

ZTGNO 

number  of  target  depths 

c 

ZTINC 

FT 

target  depth  increment 

c 

rgt 

KYd 

RANGE  TO  TARGET 

c 

RGTNO 

number  Of  target  ranges 

c 

rginc 

KYd 

target  range  increment 

c 

T  ARG 

,LT,  0.0  TARGO  IS  A  TARGET  LINE 

c 

COMPONENT  IN  DB/MICROPASCAL 

c 

.GE.  0.0  TARGO  IS  CONST  IN  TARGET 

c 

spectrum  level  equ 

c 

targo 

level  of  line  COMPONENT  in  DB  OR 

c 

const  in  SPECTRUM  equ 

c 

targi 

COEFF  OF  ALOGIO(HZ) 

c 

TARG2 

COEFF  OF  ALOG10(HZ)**2 

c 

TARG3 

COEFF  OF  ALOG10(HZ>**3 

C 

c**  OUTPUT  AND  PLOT  CONTROL  CARDS  ** 
C 


-1 

1 

1 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


OUTPT  OUTPUT  DATA  CONTROL  PARAMETER  -1 

.LE.  0  NO  VIRT  DIR  PRINT  OUT  -  ONLY 
TOTAL  SECTOR  LEVEL 

SNPLT  .LE.  0  NO  NOISE  PLOT  0 

hFLAG  .GE.  0.0  REQ  PLOT  DATA  CARDS  0 

ELIM  PLOT  DATA  THAT  00  NOT  CHANGE 
SSPLT  .LE.  0  No  VELOCITY  PROFILE  PLOT  0 


three  additional  cards  required  after  pause  for  noise  plot 
(Only  one  set  for  each  freg  and  depth  loop) 

type  format  EXAMPLE  DATA  CARD  (IN  ORDER) 

title  40ai  ambient  noise  directionalitys’  note  (*■) 

LOCATION  13a1  4020N  01715E 

Date  10a!  q7  22  71 


Two  ADDITIONAL  CARDS  REQUIRED  FOR  SOUND  SPEED  PLOT 

LOCATION  13a1  4020N  01715E 

DATE  10A1  07  22  71 
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**  DATA  CONTROL  CARDS  ** 


header  used  before  A  MESSAGE  statement 

pause  used  BEFORE  EACH  CONSECUTIVE  run 

EXCEPT  FIRST 

PRLOSS  USED  BEFORE  EACH  PROP  LOSS  ARRAY 

END-DATA  USED  AT  END  OF  DATA  OR  BEFORE  PLOT 

cards  that  are  at  end  of  data 

*****************************  begin  program  ************************> 


DIMENSION  FRE4(ll)  ,DPTH5(29)  ,FREq9(29) 

COMMON  / COMV/  ZX, ZT G , RMAX*PHID , DELPH , BNCS *■ OUTPT , V/SPD, DGREH, SNPLT , S 
lSPLT»ALAT»BERNGfHOUR.ACTlV*RAINfHFLAG»ZTGlfZTGNOrZTINC»RGT»RGTNO,R 
2GINC,TARG,TARGO,TARG1,TAKG2.TAR63,SHIPD,SHLFR,SHLFS/COMA/AVELP<30) 

3. AHB(30)  ,ARESP(30)  »ASHIP(30)  ,  ABW(30)  .APROP(30)  .ADANFOO)  ,NUMV,NUmH 
4  ,  NUMR,  NUmS  , NuMBW , NUMP*  NUmF/COMX/ CX  r  ZBM  *  PHlC  *  ALPHAC  ,  BI  ON 

COMMON  HLOS(6»60) *RNG(6*bO) »XANGC300) »YDB(300) , JPli JCNT.WlNDI ,TOTL 
H  »T0TLN»151» ITST,NTST,THRML, THRDB , CD, ZB1 , VFL AG, BANDL ( 11 > ,FRE3(U) , 
2DB(ll) ,BANDW.LOCAT(3) , DATE < 2 ) , AVEL1 ( 30 ) ,A5HI1<30) ,TDB(180) ,TDEG«18 
30 ,KPl,TRECL,TLEv *TGDP,RTG, 0(300) »W(300) »TYTLE(B) iCAPl(B) »CAP2(l3) 

4, CAP3(13> ,CAP4(l3>  ,CAP5(i3)  ,CAP6 (8 ) , CAP7 (8 ) , NMS2 , FREQ.PHINC 

DATA  IPOZ/5HPAUSe/»IPOP/5HPRLOS/ 

DATA  ZTG/20 • / * SHLFR/900 ./ , SHLFS/Q • / , RMAX/500. / » OUTPT/— 1 « / 

DATA  RADCON/57. 2957795/. SHIPD/5. 2/ rNUMS/9/ 

DATA  SNPLT »SSPLT, PHID » DELPH. eNCS,WSPD/.0. .0* .  0 , 170 • , 1 0 . *  5 . 0/ 

DATA  (AHB(I) , 1=1 , 10 ) /0 , . iO . r 20 . * 13.6, 40 . , 14. 5, 50 , r 16. 5. 90 » , 17, / 
DATA  DGREH/36 0. /,AL AT/40. /.NUMH/10/.HOUR/1500./ 

DATA  ACT  I V/.0/*B ION/ 0.0/, RAIN/0.0/ 

DATA  (ABwd)  .1  =  1, 4)/0. 0,0. 0, 1.0*0. 0/ 

DATA  NUMBW/4/,ZTGl/-1.0/,HFLAG/OtO/iZTGNO/1.0/,RGTNO/1.0/ 

DATA  (ARESP(I) , 1=1 ,4 ) /-9Q  .  , 0 . 0 , 90 . , 0 . 0/ 

DATA  NUMR/4/ 

DATA  (A5HIP(I) ,1=1,9 )/93, 1.94*154,5,215,2.92,277.2.1,337/ 


10  CALL  REAOIN  (IPOS) 

hflgi=hflag 

JDPTH=0 

JFREQ=0 

SORT  FREq  AND  DEPTH  ARRAY  INTO  Two  SEPARATE  ARRAYS 

DO  15  KDF=1*NUMF 

IF  (ADANF(KDF)  .LE.  2.)  GO  TO  H 

JDPTH=JDPTH+1 

DPTH5(JDPTH)=ADANF(KDF) 

GO  TO  15 

11  jfreq=ufreq+i 

FREQ9 ( JFREQ) =ADAnF ( KDF) 

15  CONTINUE 

DO  160  L3=l *  JDPTH 
ZX=DPTH5(L3) 

DO  160  L4=l , JFREg 

FREQ=FREQ9(Li;  !- 

IF  (L3+L4  .GT.  nFLAG=— 1. 

PTEST=0.0 
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store  for  Print  out  arrays  that  change 

IF  IVFLAG.GT.n.O)  Go  To  20 
CO  20  I  =  1.N||MV 
AVEL1(I)=AVfLP(I> 

20  CONTINUE 

IF  (L3+L4  .FQ.  ?.)  CALL  PRINTS 
JCNT=0 

totli=o.o 
151  =  0 

HOR2I=1.0F-lB 

IF  ( FREq • LT . » 0 1 « OR ■ FREQ • GT * • 5 )  GO  TO  170 

MULTIPLE  RUn  CHECK 
IF  (VFLAG.En.0.0)  CALL  ERTHc 

CALCULATE  ll  BANDWIDTH  FREliUENCIFS 
.  BANDW=AHW(Ni.|MPW-l.)  ■  ■  .  . 

FSTftT=.5*<-RAND*+SOKT(BANDW*BANDW+4.0E+6*FREO*FREQ> ) 

DO  JO  1  =  1 r 1 1 
BB5=r-l 

FREa (I) =BB5*B ANDW* . 1 
FRE3(I)=FSTrT+FrE4(a ) 

BANCL ( I ) =FUmU ( ABW  »FrE4 ( I ) *NUMBw) 

30  CONTINUE 

determine  Biological  source 

IF  (ACTIV.Gt.0.0)  CALL  BIO 

calculate  Surface  noise  source  intensities 

CALL  SURF! 

ALPHAC=ALFA(40.0.FR£0) 

H0R<iI  =  1 . 0E“1  8 

IFlsHLFs  .Eo.  0.0)  GOTO  32 

CALCULATE  I nTENSI TY/STeR ADI  AN  FOR  HOR  NOISE  DUE  TO  OIST  SOURCFS 
CALIBRATE  HoR  NOISF  WITH  VLAM  DATA  (ASSUME  10  SHIPS  OVER  shelf 
AT  hOO  NM) 

RDIsT  =  ShLFR  +  2D0(). 

SHLFD=10.*Al0G10(1.+SHLFS> 

RL0S=-86.+ls.*AL0Glu(RDlST)+ALPHAC*(RDrST*.0nl) 

00  31  J=lrl| 

BB1=AL0G10(fRF3( J)+b. >-.25 
SL1=-1 .*BEl  +  l.lfS 
SL2=3. 3*861-6.274 
C0RLF-3B. 

31  Db(  J) =-10.-4 0.*ALOG10(1u***SL1+10.**SL2>-RLOS+ShLFD+CORL 
call  BWIDTH ( HORZ I ) 

•  CONVERT  TO  INTENSITY/HORIZ  DEG 

HORZI=HORZI*DGRFH*. 0027777 
3?  CONTINUE 


H0Rd6=10.*AlOG10(HOkZI/. 109662)+fl4.8856 

PHICH  IS  THf  HALF  Angle  WIDTH  OF  SOUND  CHANNEL  SEEN  BY  RECEIVER 
SORsV  IS  SOUND  SPEFO  AT  NOISE  GENERATOR  DEPTH 
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CXrFUNU  CAVEUP » ZX ,NUMV ) 

PHICH=0.0 

S0RSV=AMINl(AVELp(2) .AVELP(NUMV) ) 

IF (S0R5V  .GE.  CX)  PHICH=aCOS (CX/sORSV) *RaDC0N 

SKIP  RAYTRACE  IF  PROP  LOSS  READ  IN 
1F(IP0S  .EQ.  IPOP)  GO  TO  175 

IF  ( SNPLT .LE  •  0  ■  0  .AND.  OljTPT  .GT.  0.0)  PRINT  180 


CS=AVELP(2> 

PHIC=0.00001 

search  For  maximum  near  surface  sound  speed 

CM=CS 

DO  40  I=lfNJMV*2 

IF  <AVELP(I) .GT. 3000.0)  60  TO  50 
C=AVELP ( I  +  l ) 

IF  (C.GE.CM)  CM=C 
40  CONTINUE 

50  IF  (CM.UE.CS)  GO  TO  60 

IF  (CX.UT.CM)  PHiC=ACOS(CX/CM)+0. 00001 
60  PHISC=RADCON*PHIc 

ZBM=AVELP(NUMV-1) 

CB=AVELP(NUMV) 

OBTAIN  LIMITING  ANGLE  FOR  BOTTOM  REFLECTIONS 
PHILIM=0.0 

IF  (CX.LT.CB)  PHiLIM=ACOsCCX/CB) 

INITIALIZE  PARAMETERS 
PHI  LOOP 

CALL  PRERAY  (0.0,0) 

IBB=0 

ICZ=0 

OBTAIN  INITIAL  ANGLE  INCREMENT 

PHST=ABS(PHID)+0.5*DELPH 

PHNO=PHST-DELPH 

PHINC=DELPH/20.0 

PHIzPHST+PHINC 

PHND-AMAXI (PHISC,PHND)/RAOCON 

70  phi=phi-phinc 

PHIS=PHI/RADCON 
IF  (PHIS.LT.PHND)  GO  TO  60 
IF  (ICZ.EQ.1)  GO  TO  90 
IF  (PHIS.GT.PHILIM)  GO  TO  90 
IBB=1 

PHIS=PHILIM+.0001 
GO  TO  90 
00  CONTINUE 

ifcicz  .eq.  h  goto  100 


DRAW  TROUGH  IN  NoISE  CURVE  NEAR  HORIZONTAL  IF  PHICH  .GT.  .5 
AND  GET  DISTANT  SOURCE  (H0R)  NOISE  COMPONENT 
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JCNT= JCNT+1 
I51=JP1+1 

XANG  ( I 51 ) =-PHICH* . 7 
XANG<I51+1)=PHICh*.7 

YOB(I51)=10.*ALOg10<10.**(THRDB/10. )+10.**<HORDB/10. ) )+100. 

I F ( PH ICH  ,LT.  .5)  rDB(I51)=Y0B(JPl) 

Y0B(I51+1)=YDB(I51) 

T0TLI=T0TLI*H0RZl*2.*PHICH 

jPl=JPl+3 

c  initialize:  cz  folding  range  check 
call  prepay  (0*0,0) 

ICZ=1 

IF(PHnD  .gt.  PHILIM)  goto  110 

PHINC=PHlLIM-PHIc 

IF  (PHINC.LE.0.0)  GO  TO  110 

PHINC=-1.43£4*PHiNC 

PHI=(PHIc+Q.001)*RADCON 

GO  TO  70 

90  CALL  RAPaTH  < ICZ, PHIs*DELR> 

IF  (ICZ.EQ.l)  GO  TO  100  - 

IF  (IBB. EG. 1)  GO  TO  60 
IF  (DELR  ,GT.  1.)  PHINC=(J.5*PHINc 
GO  TO  70 
C 

100  CONTINUE 

IF  (PHIS.LT. PHILlM)  GO  TO  70 
110  CONTINUE 

CALCULATE  OMNI  NoISE  LEVEL 

IF  (ToTLI.LE.0.0)  GO  TO  171 
T0TLN=10.0*ALOGlO(TOTLI ) +184.8856 

hordb=horob+ioo. 

FRE1=FREq*100C.O 
PRINT  200*  DGREH,TOTLN*FREl*ZX 
C ANGL-2 • *PHICH 

IF ( SHLFS  .GT.  O.o  .AND.  PHICH  .GE.  0.)  PRINT  220*HORDB*CANGL 

bypass  target  if  prop  loss  read  in 

IF(PTEST.EQ.I.O)  GO  TO  120 


IF  (ZTG1)  120*130*130 
120  IF  (SNPLT.GT.0.0)  CALL  Pi_T  (PTEST) 
GO  TO  150 

target  range-depth  loop 
130  nt=ztgno 

NR-RGTNO 

htst=hflag 

DO  140  I=1*NR 
DO  140  J=1*NT 
YT=J-1 
TR=I-1 

rtg=rgt+yr*rginc 

TGOP=ZTG1+YT*ZTInC 

CALL  TARGET 

IF  (SNPLT.GT.0.0)  CALL  PlT  (PTEST) 
HFLAG=-1.0 
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140  CONTINUE 
HFLAS=HTsT 
150  CONTINUE 
C 

c  PLOT  SOUND  SPEED  PROFILE  IF  CALLED  FOR 

IF  (SSPLT.&T.O.O )  CALL  SSPLOT 
C 

160  CONTINUE 

HFLA6=HFLGl 

IF  (IPOS.EQ.IPOZ)  GO  TO  10 
C 

STOP 

C 

170  PRINT  210 »  FREQ 
GO  TO  16o 

171  PRINT  19o 
GO  TO  150 

C  READ  PROP  LOSS  INSTEAD  OF  CALCULATING  IT 

175  CALL  AUXPR ( IPOS'  PHICH'HONZl ) 

J  PTEST=1.0 

GO  TO  1 1 0 


180  FORMAT  (//•///57X,’OUTPUT*///24Xt3HPHIf0X,4HDPHl»8X»5HMAX  R»5X»8hLO 
1SS(DB) *6X»5HDB  Up»6X,7HD8  DOWN»3x»l0Hl  STER  COR/) 

190  FORMAT  (•  TOTLI  ZERO  OR  NEGATIVE  -  GO  TO  NEXT  CASE’) 

200  FORMAT  (//6X,F7.1>'  UEG  SECTOR  LEVEL  =  *,F7.2,'  OB  F0R'»F6.1»»  HZ 
1 AT  1 rF7.1*  *  FT  DEPTH'/) 

210  FORMAT  (7H  FREQ  = > F 1 0 . 4 t 2BHKHZ  OUTSIDE  10“500  HZ  REGION) 

220  FORMAT  (bXr'SOFAR  CHANNEL  PER  STr  NOISE  LEVEL  =*F6.1«'  DB  FOR’FS.l 
1  *  36H  VERT  DEG  TOTAL  CHANNEL  'LOOK’  ANGLE/) 

end 


FUNCTION  ALFA  <T,F> 

INPUTS 

T  deg-f  temperature 

F  KHZ  FREQUENCY 

ALFA  IS  THE  ABSORPTION  COEFFICIENT 
FSQ=F*F 

fcube=f*fsq 

FLOG=ALOg(F> 

TC=5.0* (T-32.0)/9.0 

FT=(6.0*TC+118.0)/(TC+273.0) 

FT=21,9*e;XP<2.3026*FT) 

ALFAT=0.0542*EXP(1.5*FLO&) 

FCUT=FCUBE/32.76a 

IF  (FCUT.LT. 1.0>  GO  TO  10 

aFl=l,0/(l.0+FCUT) 

oF2=1.0-BF1 

GO  TO  20 

10  bF2=l. 0/(1. 0+1. 0/FCUT) 

BF1=1.0-BF2 

20  ALFAM=FSo/FT+(0.fe505317*FT/ (FT+FsQ/FT) +0,0 26847) 

ALFA=BF1+ALFAT+6f2*ALFAM 

return 

END 
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SUBROUTINE  AuXPR ( IPOS.PHxCH.HORZi ) 

THIS  SUBROUTINE  CALCULATES  AMBIENT  NOISE  WHEN  PROP  LOSS  IS  READ 
IN  AND  NOT  CALCULATED 


COMMON  /COMV/  ZX,ZTGrRMAXfPHlD»DELPH*BNCS»OUTPTrl»SPD»DGREHf SNPLT,S 
lSPLT,ALAT*BERNG*H0UR,ACTiV,RAlNfHFLAG»ZTGl»ZTGN0,ZTlNC»RG7rRGTN0»R 
2GlNC*TARG*TARG0*TARGl»TARG2f TARG3»SHIPDrsHLFRrSHLFS/COMA/AVELP(30) 

3  t AHB ( 30 ) ,ARESPf30) ,ASHIP(30) ,ABW(30) rAPRoP(30) , AD ANF ( 30 > . NUMV r NUMH 

4  r  NUMR  >  NUmS  »  numbw , NUMP  f NUmF/COMX/CX  »zbm»phic*alphac»bion 

COMMON  HLOS ( 6 » 60 ) fRNg(6*60) rXANG(300) t YDB<300) fJPl ► JCNT t WIND I r TOTL 
1I»T0TLN» 151 » ITST ,nTST  fTHRML  f  THRDBf  CD i ZB1 »  VFLAG  f B ANDL (11) rFRE3 (11) f 
2DB (11) f  BANDW  f  LOC AT (3 ) fDATE(2)  f A VeLI < 30 ) f ASHI 1 < 30 ) , TDB ( 180 >  f TDEG ( 18 
30) fKPI fTRECLfTLEVfTGDPfRTGfQ<300) fW(300) r  T YTLE ( 8 ) rCAPl (8) »CAP2( 13) 
4 i C AP3 ( 13 ) r  C AP4 ( 1 3 ) »CAP5<13) fCAP6(8>  fCAP7(6) f NMS2 t FREQ  *  PH  INC 


DATA  IPOZ/ 'PRLOS i / 


INITIALIZED  PARAMETERS  ARE  THE  FOLLOWING 
T0TLI=1.0E-18 

T0TLI=T0TLI+H0RZl*2.»PHICH 
UP1  =  0 

HORDB  =  10.*AL0Gi0(HoRZl)+64.8856 
IF  (PHICH  .LT.  .5)  GO  TO  1 
Q(1)=-PHICH*.7 
Q(2)=-Q(l) 

W(1)=10.*ALOG10 (10.**(THRDB/10. )+10.**(HORDB/10. ) )+100. 
W(2)=WU) 

XANG ( 1 ) =1 . 

XANG(2)=1. 

JP1=2 

1  CALL  READIN(IPOS) 

CALCULATE  WIND  NoISE  DIRECTIONALITY  FACTOR  DIRS 

CX=FUNU(AVELPrZX,NUMV) 

CS=FUNU(AVELP»ZTg»NUMV) 

COSPH=  CoS(.01738*APR0P(l))*CS/CX 
IF(COSPH.GE.l. )CoSPH=l. 

PHI=ACOS(COSPH) 

EXPN=5.729#PHI 
IF (PHI ,LT. .3491 )EXPN=2.0 
DCOEF=2.*FREa-.02 
IF(FREQ.GT..5)DCOEF=1.0 
DIRS=1.  +  (SIN  LPHI ) ++EXPN-1 . ) * ( DCOeF ) 

C 

RI-APR0P(3) 

SHI=FUNS(ASHU,RifNMS2) 

PRI=APR0P(4) 

RECTI=1.0E-1B 
I=-l 
4  1=1+1 
NS=2*1+1 

IF (ASH II (NS) .LE*APR0P(3) )  GO  TO  4 
NP=5 
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c  select  range  integration  interval 

10  IF ( APROP ( MP ) -ASH 1 1  { NS ) ) 30 *20, 40 
20  RTEST =1  •  0 
NS=NS+2 

30  HIP1=APRoP(NP) 

PRIP1=APR0P(NP+1) 

SHIP1=FUNS(ASHU,RIP1*NMS2) 

NP=MP+2 
GO  TO  50 

40  RIP1=ASHH  (NS) 

PRIP1=FUNU( APROP. RIP1 ,NUmP> 

SHIP1=ASHI1(NS+1) 

NS=NS+2 

lF(RTEST.EQ.l.O)  NP=NP+2 
50  CONTINUE 

CALCULATE  TOTAL  SURFACE  NOISE  intensity  density  FUNCTION 
dni=a*range+b 

kdif=ripi-ri 

ifirdif.eo.o.o)  go  to  ao 

SDIF=SHIP1-SHI 
A=SDIF/RDIP 
bwind=windi*dirs 

bSHlP=SHI-RI*A 
b=BWlND+BSHIP 

dni=shi+bwind 

DNIP1=SHIP1+BWINd 

calculate  prop  loss  intensity  reduction  factor 

FROM  PROP  LOSS  =  E*RaNGF  +  F  ,  WHERE  F  INCLUDES  G  The  LOSS 
DUE  TO  THE  BEAM  PATTERN 

PDIF=PRIP1-PRI 
E=PDIF/RDIF 

G=FUNU<ARESP»APRoP<1) *NUmP) 

F=PRI-RI*E-G 

IF ( ABS ( E )  .GE.  .0001)  GO  TO  60 

RECI=6.282*lQ.**<-. 1*F>*( . 3333*A* (R IP1**3-R I **3 ) +.5*B* ( RIP1*R IP1 
1-RI*RD  ) 

GO  TO  70 

60  RECPE=10./(E*2.3o25) 

RECI=-6.282*RECPe*<10.**(-.1*PRIP1)*  DNIPl* ( R1P1+RECPE ) -10 .** 
1(-.1*PRI)*<DNI*(RI+RECPE) ) ) 

IF(ABSU)  .GT,  .0001)  RECI=RECI-6.282*A*RECPE*RECPE*(10. *«■(-. 1* 
1PR1P1)*(RIP1+2.*RECPE)-10.**(-«1*PRI>*‘RI+2.*RECPe) ) 

CALCULATE  LEVEL  RECEIVED  FROM  HORIZ  SECTOR  AND  TEST  FOR  LAST  RANGE 

70  RECTI =RECTI+.0175*DGREH*ABS (RECI ) + (THRML+B ION) * APROP (2) 

80  RI-RIPl 


SHI=SHIP1 

PRI=PRIP1 

IF(NP.LT.NUMP)  Go  TO  10 
C 

C  STORE  LEVEL  AND  ANGLE  FOR  OUTPUT  AND  GET  TOTAL  OMNI  INTENSITY 
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JP1=JP1+1 

W( JP1)=10.0*ALOG10<RECTI/APROP<2) )+lB4.8856 
Q(JP1)=APR0P(1) 

XANG(JP1)=APR0P(2) 

TOTLI=TOTLI+RECTl 
IF(IP0S  .EQ.  IPOZ)  GOTO  1 

SORT  AND  GET  MINIMUM  LEVEL 

UP1M1=JP1-1 
100  DO  110  K=lrJPlMl 
L=K+1 

IF<q(|_) .GT.Q(K) >G0  TO  110 
DUMYQ=Q(L) 

QUMYW=W(L) 

OUMYXrXANG(L) 

&<L)=Q(K) 

W<L)=W(K) 

XANG<L)=XANG(K) 

G)<K)=DUMYQ 
ml  ( K )  =DUMYW 
XANG ( K ) -DUMYX 
110  CONTINUE 

JP1M1=JP1M1-1 

IF ( JP1M1 .  GE ,2 ) GO  TO  100 

PRINT  OUT  IF  NO  PLOT  REQUIRED 

IF(SNPLT.GT.0.0)gO  TO  130 
PRINT  200 
DO  120  K=1,JP1 
120' PRINT  2 1 0  ’  Q  ( K  ) »W(K) 

GO  TO  150 

130  DO  140  K=1,JP1 

140  W<K)=W(K)-10,*ALoG1O(XANg(K) ) 

150  CONTINUE 

200  FORMAT(//llX, »PHi* »4X* ‘LEVEL  ( DB ) • ) 

210  FORMAT (Fl5»l,F10,l) 

RETURN 

END 


subroutine  bio 

this  SUBROUTINE  COMPUTES  THE  CONTRIBUTION  OF  BIOLOGICAL  NOISE 

COMMON  /COMV/  ZX,ZTG,RMAX»PHID*DELPH,BNCS'OUTPT,WSPD,DGREH,SNPLT,S 
1SPLT, ALAT,BERNG* HOUR, ACTiV, RAIN, HFLAG»ZTg1.ZTGNO,ZTINC#RgT,RGTNO,R 
2GlNC*TARGrTARGO*TARGl»TARG2,TARG3,SHIPD,sHLFR,SHLFS/COMA/AVELP(30) 

3, AHB<30) ,ARESP(3o) ,ASHIP(30) ,ABW(30) ,APROP<30) ,ADANF(30) ,NUMV»NUmH 

4 ,  NUMR , NUmS , NuMBW , NUMP , NUmF/COMX/CX  » ZBM »  PH I C  »  ALPHAC , B ION 

COMMON  HLOS«6*60) ,RNG(6»60) ,XANG(300) ,YDBl300) ,JPl*JCNT»WlNDI,TOTL 


1I'T0TLN»I51*ITST,NTST  r  THRML , THRDB , CD , ZB1 , VFLAG  *  BANDL (11)  iFRE3(11)  , 
20B(11) ,BaNDW,L0CaT(3) , DATE (2) .AVELl (30) »ASHI1 (30) ,TDB(180) ,TDEG< IB 
30) , KPl , TRECL , TLEV,TGDPrRTG*Q(300) ,W<300) ,TYTLE(B) ,CAP1 (3) rCAP2(13) 
4,CAP3(13J ,CAP4(l3) *CAP5<13) ,CAP6(8) , CAP7 (8 ) , NMS2 , FREQ* PHINC 


50 


ooo  ooo  o  o  o  o  o  o  o  noon  ooo 


ACT=ACTlV/5.0 
DO  10  J=l.ll 
BB1=ALOG10(FRE3(J) ) 

8  IS  THE  DAILY  FLUCTUATION 

B=. 00175* ( 100. -FrE3(j) >*SIN( . 002618* ( HOUR-300. ) ) 
Bl01=-26.6*BBl+2o. 

BI02=28.*BBl-52.64 

10  QB( J)=-35.4-10.*aLOG10< 1o.**( . 1*8101 ) +10 .**(. 1*BI02) )+B*ACT 
l+2.*ACTIv 
CALL  BWIOTH  (B2> 

CONVERT  TO  INTENsITY/VERT  DEG/H0R  DEG 

BION=ACT*B2*DGREH*l .5432E-5 

RETURN 

END 


SUBROUTINE  BwIDTh  (DINT) 

THIS  SUBROUTINE  INTEGRATES  THE  SIGNAL  INTENSITY  OVER  BANDWIDTH 

ADD  BANDWIDTH  RESPONSE  To  NOISE  SIGNAL 

REAL  LEVBW ( 1 1 1 i Ml f MIP1 

COMMON  /COM V/  ZXt2TG»RMAXrPHlD»DELPHiBNCS»OUTPT(wSPD»DGREH»SNPLT,S 
lSPLT.ALATfBERNG»HOUR»ACTlV.RAlNfHFLAG*ZTGlrZTGNO,ZTlNC*RGT»RGTNOrR 
2GlNC»TARGfTARGOrTARGl»TARG2»TARG3,SHIPD»SHLFRrSHLFS/COMA/AVELP(30) 
3 1 AHB (30) »ARESP(3o) . AsHIP (30 ) . ABW ( 30) . APROP ( 30 ) » ADANF ( 30 ) » NUMV » NUMH 
4 » NUMR ► NUmS . NUMBW , NUMP r NUmF 

COMMON  HL0S(6#60)  » RNG  ( 6 »  60  )  .  XANG  ( 300  )  .YDBO00)  *  JPl  .  JCNT  »  WlNDI  r  TOTL 
1 1 f TOTLN . 1 51 . ITST , NTST  r  THRML. THRDB  .CD.ZB1* VFLAG .  BAnDL (11) .FRE3 ( 1 1 ) , 
2DB (11) f BANDW  *L0CaT (3) *  DATE (2) » AVEL1 ( 30) iASH 11(30) .TDB( l8o) *  TDEG ( 18 
30) »KPlfTRECL»TLEV»TGDP»RTG»Q(300) ,W(300) .TYTLE(8) .CAP1<8) .CAP2U3) 
4,CAP3(13) »CAP4(13) .CAP5<l3) rCAP6(8) »CAP7(8) .NMS2.FREQ.PHINC 


DINT=0.G 

OBTAIN  effective  LEVELS 
DO  10  1=1.11 

10  LEVBW(I)=DB(I)-BaNDL(I) 

INTEGRATE  AND  CONVERT  TO  INTENSITIES 

DO  30  1=1.10 
FI=FRE3( I ) 

FIP1=FRE3(I+1) 

Ml  =  (LEVBw  < I  +  l )-LeVBW ( I ) ) / (10 .0*AlOG10 (FIP1/FI ) ) 
BI=LEVBW(I)-(ALOG10(FI)/ALOG10(FIP1/FI) )*(LEVBW(I+1)-LEVBW(I) ) 
SMB=3,2467E~9*10.**(BI/10. ) 

IF  (MI.EG.-1.0)  GO  TO  20 
MIP1=MI+1 . 0 

DIN= ( SMB/MIP1 ) * (FIP1**MIP1-FI**MIP1) 

GO  TO  30 

20  DIN=SMB*ALOG(FIPl/FI ) 

30  DINT=DINT+DIN 
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c  normalize  intensity  to  ONE  HZ 

DlNT=DINT/BANDW 

C 

RETURN 
END 

SUBROUTINE  erthc 

THIS  SUBROUTINE  CORRECTS  for  earth  CURVATURE  EFFECTS 

common  /CCMV/  ZX,ZTG.RMAX*PHID»DeLPH.BNCS>OUTPT»wSPD»DGReH#SNPLT>S 
1SPLT  > ALAT r  BERNG  »H0UR  *  ACTl V » RAIN#  HFLAG  f ZTgI » ZTGNO#  ZTI NC  »  RgT » RGTNO f  r 
2GlNC»TARG»TARGO»TARGI»TARG2.TARG3fSHIPD»SHLFR,SHLFS/COMA/AVELP(30) 
3* AHB<30) ,ARESP(3o) » ASH IP ( 30 ) » ABW ( 30 ) »APROP(30) » ADANF < 30 > » NUMV > NUmH 
4 » NUMR » NUmS » NUMBW , NUMP » NUmF/COMX/CX * ZBM » PH IC » ALPHAC »B I ON 
COMMON  HLOS(6»60) »RNG<6rbO) >XANG(300) rYDB<300) r JP1» JCNTrWlNDI »TOTL 
H  »TOTLN# I51» lTST,NTSTrTHRML»THRDB»CD.ZBl,VFLAG*BANDL( II) rFRE3(ll) , 
2DB ( 11 ) ,BANDW,LOCAT(3) .DATE (2) ► AVEL1 ( 30 ) t ASHI1 <30)»TDB(18Q)»  TDEG (18 
30) »KPl»TRECL*TLEv*TGOP*RTG»Q(30Q) r  W ( 30  0 ) t TYTLE  tfl) »CAPl(0)f  CAPS (13) 
4  > C AP3 (13) t CAP4 (13) fCAP5(l3),CAP6(B) »CAP7(8) rNMS2f FREQ*PHINC 

DATA  RADCON/57. 2957795/ 

DATA  AB/43642464.29E7/»ASQ/43789029.Q0E7/>BSQ/43496390.14E7/ 

STORE  ORIGINAL  BOTTOM  DEPTH  FOR  PLOT  DATA 
ZB1=AVELP(NUMV-1) 

earth  CURVATURE  CORRECTION 

JPAIRS=NuMV/2 
alat=alat/radcon 

SINL=SIN (ALAT) 

COSL=COS(ALAT) 

sinlq=sinl*sinl 

COSLQ=COSL»COSL 

ROL=AB/SqRT(ASQ*SINLO+BSq*COSLO) 

DO  10  I=1»JPAIRS 
12=2*1 

AVELP ( 12  J SAVELP ( 1 2 ) *ROL/ ( ROL-AVELP ( 1 2-1 ) ) 

AVELP(I2-1 )=AVELp( 12-1 )*(1.+(AVELP{ 12-1)/ (2.*R0L) )+( < AVELP ( 12-1 ) /R 
10L)**2. )/3») 

10  CONTINUE 
VFLAG=1.0 
ALAT=ALAT*RADCON 
RETURN 
END 

FUNCTION  FUNS(A*X>N) 

this  function  interpolates  the  shipping  histogram 

DIMENSION  A ( 1 ) 

ILAST=N 
1=1 

. IF(X.GT.A(I).AND.X.LE.A(ILAST))  GO  TO  10 
FUNS=0.0 
RETURN 
10  1=1+2 

IFIX.LE.A(I)  .AND.X.GT.AQ-2)  )  GO  TO  20 
GO  TO  10 
20  FUNS=A(I-1) 

RETURN 
END 
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FUNCTION  FUNU  <A,X#N) 

DIMENSION  Ail) 

ILA5T=N-1 

1=1 

IF  (X.GT.A(I) )  GO  TO  10 
FUNU=A<2) 

RETURN 
10  1=1+2 

IF  (X.LE.A(I))  Go  TO  20 
IF  (I.LT.ILAST)  GO  TO  10 
GO  TO  30 

20  IF  <A(I).EQ.A(I-2)>  GO  TO  30 
IF  (X.EQ.A(I))  Go  TO  30 

FUNU=A ( I -1 )  +  ( A  U  +1 ) -A  < I -1 ) ) / ( A (I ) -A (1-2 ) ) * { X-A ( I -2 ) ) 
RETURN 

30  FUNU=A<I+1> 

RETURN 

END 

SUBROUTINE  HLOSS  (PHIE.R.S.DRDP>GRAD»HkSR) 

THIS  SUBROUTINE  COMPUTES  ABSORPTlON(HK) .  SPREADING (HS) . 
AND  REFRACTION (Hr)  LOSSES. 


COMMON  /COMX/  CX , ZBM.PHIC . ALPHAC «  BION/COmCT/CV»  T ANPHS 

DATA  FMIN/1.0E-6/fFMAX/1.0E12/ 

IF  (ABS(PHIE) .LT.0.0001)  GO  TO  10 
FOCUS=ABS(CV/CX*r*SIN(PHiE)*DRDP) 


GO  TO  20 

10  FOCUS=ABS(CV/CX*R*TANPHS*CV/GRAD/3000.0) 

20  IF  (FOCUS. GT4FMlN. AND. Focus. LT.FmAX)  GO  TO  30 
PRINT  40.  PHIE.R.DROP.GRAD 
CALL  EXIT 

30  HKSR=60. 0+1 O.O+AlOGIO (FOCUS )+ALPHAC*S 
RETURN 

40  FORMAT  (5X. 'PHIE.R.DRDP.GRAD' .4E12.6) 

END 


FUNCTION  HSURF  (n.NP) 

THIS  SUBROUTINE  COMPUTES  SURFACE  REFLECTION  LOSS 

N  =NUMBER  OF  reflections 
np=path  type 


NP 

1 

2 

3 

4 

5 

6 


PATH 

DU 

DD 

UU 

UD 

DP 

SR 


REFL'S 

N-l 

N 

N 

N+l 

0 

1 
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C 


DATA  HS/0.5/ 


GO  TO  (10, 20,20.30. 4o»50) .  NP 

10  ens=n-i 

GO  TO  60 
20  ens=n 
GO  TO  60 
30  ENS=N+1 
GO  TO  60 
40  ENS=0.0 
GO  TO  60 
50  ENS=1.0 
C 

60  HSURF=ENS*HS 
C 

RETURN 
END 

SUBROUTINE  PLT  (pTEST) 

THIS  SUBROUTINE  PLOTS  THE  DIRECTIONAL  AMBIENT  NOISE 

COMMON  /COMV/  ZX.ZTG.RMAX. PH ID,  DELPH .BNCS »  OUTPT . WSPD .  DGREH , SNPLT  . S 
ISPLT.ALAT.BERNG .HOUR, ACTIV. RAIN .HFLAG.ZTGl.ZTGNO.ZTlNC.RGT.RGTNO.R 
2GINC.TARg*TARG0.TARG1.TarG2.TARG3.SHIPD.SHLFR.SHLFS/C0MX/CX,ZBM.PH 
3IC.ALPHAC.BI0N 

COMMON  HLOS<6»60) ,RNG<6,bO) .XANG(3Q0>  »YDB<300) . JPl. JCNT.WINDI.TOTL 
1I.T0TLN. I51»ITST,NTST. THRML. THRDB .CD . ZB1 . VFLAG. BANDL (11) »FRE3(ll) . 
2DB(11) .BANDW.L0CAT(3).DATE(2) .AVELK30) »ASHI1(30) ,TDB(180) .TDEGC18 
30)  .KPI.TRECL.TLEV.TGDP.RTG.QOOO)  »W(300>  ,TYTLE(8)  ,CAPK8)  .CAP2C13) 
4 . CAP3 ( 13) »CAP4(l3)»CAP5(l3) ,CAP6(8> .CAP7 18) .NMS2. FREQ. PHI NC 


integer  PN 

DATA  BERNG/270./.PN/0/ 

CT=CD+80.0 

IBRrBERNG 

idgr=dgreh 

DBMAX=-l000. 

IF  <hflag.lt.o.o  .OR.HFLaG.GT.O.O  .and.numf.gt.d  GO  TO  10 
READ  (5.160)  (TYTLE(I) .1=1.8) 

READ  (5.160)  (LOCAT(I) .1=1.3) 

READ  (5.160)  (DATE(I) .1=1.2) 

10  CONTINUE 

SET  UP  AnD  PRINT  OUT  HEADINGS  IF  REQUESTED 

ENCODE  (40.160. CaPD  (TYTLE(I)»  1=1.8) 

WSP1=WSPD 

FRG=FREQ*1000.0 

ENCODE (64. 220. CAP2) (LOCAT ( I ) , 1=1 . 3 ) .FRQ 
ENCODE <64. 230 .CAP3) (OATE ( I ) . 1=1 . 2 ) .ZX 
ENCODE ( 64.240 .CAp4)  wSPl.TOTLN 
ENCODE(64.250»CAP5)  ZBl.lDGR.IBR 
IF  (OUTPT  .L£.  0.)  GO  TO  15 
PRINT  20o»  (CAP2(I) .1=1.13) 

PRINT  200.  (CAP3(I) .1=1.13) 

PRINT  200.  (CAP4( I) ,1=1,13) 
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PRINT  200#  <CAP5ll> »I=1»15> 

PRINT  170 
PRINT  180 
15  CONTINUE 

IF  (PN.Le.O)  call  BALLPT  ('PEN  2  BLACKS') 
PN=PN+1 


SKIP  SORT  IF  PROP  LOSS  READ  IN 
IF(PTEST.EQ.I.O)  GO  TO  151 


SORT  DATA  AND  SET  MAX  AND  MIN  GRAPH  LIMITS 
20  N=0 

J=I51-1 

DO  30  I=1#J»2 

N=N+1 

Q<N)=XANG(I) 

IF  (YDB<I).LE.CD>  YDB«I)=CD 
30  W(N)=YDB(I) 

M-0 

DO  40  1  =  151 #  JP1 #  2 
M=M+1 

N=JCNT-M+1 

QlN)=XANG(D 

IF  (YDB(I) .LE.CD)  YDB(I>=CD 
IF  (YDB(I) .GE.CT)  YDB(I)=CT 
40  w(N)=YDB(I> 

N= JCNT 
L=I51+1 

DO  50  I=L  t JP1 #2 
N=N+1 

Q(N)=XANg(I) 

IF  (YDB(I) .LE.CD)  YDB(I)=CD 


IF  (YDBlI) .GE.CT)  YDB(I)=CT 
50  W (N ) =YDB ( I ) 

60  M=0 

DO  70  1=2  1 151 #2 
M=M+1 
N=JP1+1-M 
Q(N)=XANG(I) 

IF  (YDB(H .LE.CD)  YDB(I>=CD 
70  w  IN)=YD8  1 1 ) 

PRINT  ANGLE  AND  LEVEL  WHEN  PROP  LOSS  NOT  READ  IN  IF  REQUESTED 
IF  (OUTPT  .GT.  0.)  PRINT  190»Q(1)#W(1) # Q ( JP1 ) #W ( JPl ) 

1  =  1 

80  SUM=0 . 0 
90  1=1+1 

IF  ( I .GE. JCNT )  Go  TO  110 
SUM=SUM+ABSlQ( I)-Q(I-I) ) 

IF  (SUM-1.0)  90»100#100 


100 

L= 

JP1-I+1 

IF 

(OUTPT  .GT. 

0.)  PRINT 

190  #  Q ( I ) #W(I)iQ(L)#W(L) 

GO 

TO  80 

110 

IF 

(OUTPT  .GT. 

0.)  PRINT 

210 

CALL  BGNPL  (1) 

CALL  PHYSOR  12.0,1.5) 

CALL  TITLE  <CAP1,-100#0#C#0#0#7.o#5.5> 
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CALL  ENDgR  (1) 

CALL  PHYSOR  (2.0,1.5) 

CALL  TITLE  (1H  »-l VERTICAL  RECEIVED  ANGLE  (( >DEG <) )S' » 100 , 'NOISE 
1  LEVEL  (ODB  RE  <M)PA  /STR  HZ( ) ) S ' » 100 > 7.0 ,4.8) 

CALL  HEIGHT  (.1) 

CALL  ANGLE  (90.0) 

CALL  MESsAG ('2$*, 100, -.7, 3.2) 

CALL  RESET  ('ANGLE') 

CALL  RESET  ('HEIGHT') 

CALL  YTICKS  (2) 

CALL  XTICKS  (3) 

CALL  FRAME 
CALL  YAXANG  (0.0) 

CALL  INTAXS 

CALL  MIXALF  CL/cGREEK*) 

CALL  GRAF  (-90. 0 ,30 . 0 * 90. 0 ,CD 1 10 . 0 , CT ) 

CALL  RESET  ( 'MIXaLF') 

CALL  HEADIN  (CAPg . 100 t 1 r 4 ) 

CALL  HEADIN  (CAPa. 100. 1 .4) 

CALL  HEADIN  ( CAP4 1 100» 1  * 4 > 

CALL  HEADIN  (CAP5»100*1.4> 

CALL  CURVE  (Q»WrjPl»Q) 

PLOT  TARGET  LEVELS  IF  CALLED  FOR  AND  PROP  LOSS  NOT  READ  IN 

IF(TGDP.LE.O.O.Or.PTeST.EQ.I.O)  GO  TO  130 

KP3=3*KPl 

00  120  I=3»KP3»3 

J=I/3 

Q(I-2)=TDEG(J> 

Q(I-1)=TDEG(J) 

Q(I)=TDEG(J) 

W(I-2)=CD 

tflI-l)=TDB(J) 

W(I)=CD 


IF  (DBMAX.LT.TDB(J))  OBMAX=TDB(J) 

IF  (W(I-l) .LE.CD)  W(I-1)=CD 
120  IF  (HI(I-l).GE.CT)  W(I-1)=CT 
IZTG1=TG0P 

ENCODE (64 1 260*  CAp5)  IZTG1 »RTg 
ENCODE (64,270#CAP4)  TRECLrTLEV 
ENCODE (40  »280»  CAP6)  DBMAX 
ENCODE (40 .29o»CAp7) 

CALL  HEIGHT  (.1) 

CALL  MESSAG  (CAP7 r 100 » 4.7 r 4 .45) 

CALL  HEIGHT  (.07) 

CALL  MESSAG  (CAPS  1 100 » 4 . 3 » 4 , 25) 

CALL  MESSAG  (CAP4» 100»4. 4»4,1) 

CALL  MESSAG  (CAP6>100’4.4’3,95) 

CALL  DASH 

CALL  CURVE  (Q»W»kP3,0) 

CALL  RESET  CDASh') 

CALL  RESET  ('HEIGHT') 

130  CALL  ENDPL  (0) 

C 

RETURN 

C 

C  PRINT  OUT  ANGLE  AND  LEVEL  WHEN  PROP  LOSS  READ  IN  IF  REQUESTED 
C 

151  1=0 
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C 


152  1=1+1 
J=JP1+1-I 

IF  (OUTPT  ,GT.  0.)  PRINT  190#  Q ( J J # W ( J) #Q ( I > » W ( I ) 
lF(J-I-2>  154»153»152 

153  IF  (OUTPT  .GT.  0.)  PRINT  190#  Q ( J-l ) # W ( J-l ) 

154  CONTINUE 
GO  TO  110 


160  FORMAT  (13A5) 

170  FORMAT  (///12X . 9HD/E  ANGLE»5X#HHN0ISE  LEVeL#9X»9HD/E  ANGLE#5XllHN 
10ISE  LEVEL) 

180  FORMAT  (14X#5H(DeG) »10X#7H(DB>  * » 12X # 5H (DEG) # 10X »7H (DB)  */) 

(3X#F16.1,F14.1#F£0.1#F14.1) 

(10X# 12A5, A2) 

</7X»»  *  DB  RE  MICR0PASCAL**2/STERADIAN  HZ*/1H1) 

( • LOCATION* #4X#3A5#10X#’ FREQ  (HZ) •  »F16.1»  •*• ) 

( 'DATE* »13X#2A5»10X# »REC  DEPTH  (FT ) ♦ r Fll . 1 # • I  * ) 

(•WIND  SPEED  (KT)  •# F9. 1 » 12X #• SECTOR  LEVEL  (DB)*»F8.1»" 


190 

200 

210 

220 

230 

240 

250 


format 

FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
IS*) 

>60  FORMAT 
FORMAT 
FORMAT 
FORMAT 
END 


270 

280 

290 


»S*> 

(•BOTTOM  DEPTH  ( FT ) * »F8. 1 » 12X # *HOR  SEC  <DEG)*«I5#»  AT*#I4,» 

(• - ' # 14# •  FT  TARGET  AT • » F7. 1 # ’  KYD  S*> 

(F7.1#'  Db  REC  FR0M'»F7.1,*  DB  SOURCE*') 

(F7 » 1 » *  DB  MAX  ONE  DEGREE  AVERAGES') 

(•  TARGET  DATA  S*) 


SUBROUTINE  PRINTS 

THIS  SUBROUTINE  PRINTS  THE  INITIAL  PARAMETERS 
DIMENSION  D0UTI24) 


COMMON  /COMV/  ZX, ZTG,RMAX»PHID#DELPH#BNCS» OUTPT rWSPD»DGREH#SNPLT,S 
lSPLTrALAT»BERNG*HOUR»ACTiV#RAlN#HFLAG#ZTGl#ZTGNO#ZTINC»RGT#RGTNO,R 
2GINC#TARG»TARG0#TARG1»TaRG2»TARG3#SHIPD#SHLFR#SHLFS/COMA/AVELP<30) 
3  #  AHB ( 30 ) » ARESP ( 3o ) #  ASH IP (30 ) #ABW(30)#  APROP ( 30 ) #  ADANF ( 30 ) #NUMV#NUMH 
4 » NUMR  #  NUMS » NUMBW , NUMP  # NUmF 

COMMON  HLOS(6#60) #RNG(6»60).XANG(300) #YDB(300) . JPI#UCNT»WlNOl#TOTL 
1 1 » TOTLN #151#  ITST ,  NTST  #  THRML  »  THROB  #CD.  ZB1 ,  VFLAG  #  BANDL ( 11 )  # FRE3 ( ID# 
20B ( 11 ) >  BANDW  > LOCaT ( 3 ) #DATE<2) #AVELl(30) # ASHll ( 30 ) . TQB ( 180 ) # TDEG (18 
30) #KPl#TRECL, TLEV»TGDP»RTG. 0(300) #W(300) ,TYTLE(8) ,CAP1(8) #CAP2(13) 
4,CAP3(13) ,CAP4(l3)»CAP5(l3).CAP6(a) #CAP7(8) ,NMS2»FREO»PHINC 

DATA  BLANK/5H 


SORT  SHIPPING  HISTOGRAM  ARRAY 


I=-l 
K=-l 
1  K=K+2 
1  =  1+2 

IF(K.GE.NUMS)  GO  TO  4 

Q(I)=ASHIP(K) 

Q(I+1)=ASHIP(K+1) 

W ( I ) =20000  • 

W ( 1+1 ) =40000  • 
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2  "IF  ( ASHIP  (K+2 )  •  GE ,  0 . 0  )  GO  TO  1 

lF(ASHlP(K+2).LE.-50.)  Go  TO  3 
w(I)=-ASHlP(K+2> 

K=K+1 
GO  TO  2 

3  w(I+1)=-ASHIP(K+2) 

K=K  +  1 

GO  TO  2 

4  Q(I)=ASHIP(K) 

W(I)=20000. 

W(I+1)=40000. 

NUMQ=I+1 

max=maxo(numv.numR.numh.numq,numbw> 

PRINT  INPUT  VARIABLES 
PRINT  70 

PRINT  80,  FReQ.ZX.ZT6.PHID.DELPH,BNCS»WSPD.DGREH,bERNG.H0UR»ACTIV* 
lRAIN#ALAT»SNPLTrsSPLT*HFLAG#ZT6l.ZTGNO»ZTlNC»RGT,RGTNO»RGlNC»TARG* 
2TARG0»TARG1»TARG2»TARG3.SHIPD.SHlFR.SHLF5 


PRINT  90 
PRINT  100 

PRINT  ARRAYS 

DO  60  J=1.MAX»2 
DO  10  1=1  »  24 
10  DOUT(I)=BLANK 

IF  (J.GT.NUMV)  Go  TO  20 

ENCODE 1 20 .110* DOuT (2) )  AVELl(J).  AVELKJ+l) 
20  IF  (J.GT.NUMH)  Go  TO  30 

£NCODE<20» 110»D0uT<6) )  AHB(J) »AHb(J+1) 

30  IF  (J.GT.NUMR)  Go  TO  40 


ENCODE (20 .110. DOuT ( lo ) )  ARESP(J) .ARESP(J+l) 

40  IF  (J.GT.NUMQ)  Go  TO  50 
IF  (J.GE.NUMQ-1)  GO  TO  44 
SUM=W(J)+W< J+l) 

IF  (SUM. LT. 55000.)  GO  TO  41 
ENCODE (25. 130 .DOuT (14) )  q(J) ,Q(J+D 
GO  TO  50 

41  IF  (SUM. LT. 35000.)  GO  TO  42 

ENCODE (25. 130 .DOuT (14) >  G ( J) . Q ( J+l ) . W ( J) 

GO  TO  50 

42  IF  (SUM.LT. 15000. )  GO  TO  43 

ENCODE (25. 140. DOuT (14) )  a ( J) . Q ( J+l ) . W ( J+l ) 

GO  TO  50 

43  ENCODE (25.130 .DOuT ( lu ) )  Q ( J) . Q (U+l ) . W ( J) . W ( J+l) 
GO  TO  50 

44  ENCODE (2o » 130  » DOuT ( 14) )  Q(J) 

50  IF  (J.GT.NUMBW)  GO  TO  60 

ENCODE (20 . 110 ► DOuT (19) )  ABw( J) . ABW(J+1) 

60  PRINT  120.  (DOUT(K) .K=1.24) 

RETURN 
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70  FORMAT  (////48X.»**  INITIAL  PARAMETERS  **'//) 

80  FORMAT  (BX'FREQ  =  *  F6.3*  KHZ*5X.ZX  ='F6.1*  FT'6X'ZTS  =*  F6.1 
1*  FT»6X*PHID  =*F6.1*  DEG*5X*DELPH  =*F6.1*  DEG*/8X*BNCS  S'F6.1*9X 
2*WSPD  =*F6.1*  Kj  *5X*DGREH  =*F6,lt  DEG’5X'BERNG  ='F6.1*  DEG*5X’H0 
3UR  =*F7»2/8X*ACtIV  =  *F6. 1 . 9X*RAIN  =*F6.1*  IN/HR*  3X * ALAT  =»F6*1* 
4  DEG*5X*SNPLT  = *F6. 1 , 9X • SSPLT  ='F6.1/8X*HFLAG  =* F6. 1 . 9X *ZTG1  =  *F6 
5.1*  FT*6X’ZTGN0  =*F6.1.9x'ZTINC  =*F6.1*  FT'6X*RGT  =»F6.1*  KYD  */ 
6BX ' RGTNO  =*F6.1*9X*RGINC  =*F6.1*  KYD  *4X*TARG  =*F6.i» 9X*TARG0  =*F 
76. 1 .9X  *  TARG1  =*Fe.l/8X*TARG2  =*F6. 1 .9X* TAR63  =*F6.1»9X*SHIPD  =  *F6. 
81»9X*SHLFR  =*F6.1*  NMl'5X'SHLFS  ='F6.1///> 


1PING  HISTOGRAM*  »llX» •BANDWIDTH*/) 

100  FORMAT  (  9X»»  Ft  FT/sEC »» 10X, • DEG  OB*.12X.*OEG 
II  SHIPS  KT  FT’.sX.'HZ  DB  DOWN*/) 

110  FORMAT  (F10.1.F8.1) 

(24A5) 

( F7 • 1 , 2F6 , 1  *  F6 • 0 ) 

(F7.1.F6.1.6X.F6.0) 


DB*  *7X» '  NM 


120 

130 

140 


FORMAT 

FORMAT 

FORMAT 

END 


SUBROUTINE  RAPATh  (ICZ.PHIS.DELR) 


THIS  SUBROUTINE  GENERATES  RAY  PATHS  AND  CORRESPONDING 

propagation  losses  for  bb  and  cz  modes,  ray  path 

COMBINATIONS  ARE  OBTAINED  FOR  UP  TO  NMAX  BOTTOM 
REFLECTI0NS(0R. ReFRACTINONS)  FOR  EACH  0/e  ANGLE (PHIS). 

DIMENSION  SX  (4 )  *  ST(4> 

COMMON  /COMV/  ZX , ZTG » RMAX » PHID . DeLPH.BNCS * OUTPT . WSPD . DGREH . SNPLT , S 

lSPLT»ALAT.BERNG.HOUR»ACTlV.RAlN»HFLAGrZTGl»ZTGNO.ZTINC»RGT.RGTNO.R 
2GlNC.TARGrTAR60»TAR61.TARG2,TARG3.SHlPD»SHLFR.SHLFS/COMA/AVELP<30) 
3»AHB(30) »ARESP(3o>  »ASHIP(30) .ABW(30).APRoP(30)f ADANF(30) .NUMV.NUmH 
4,NUMR.NUmS*NUMBW,NUMP.NUmF/C0MX/CX.ZBM»PHIC»ALPHAC.BI0N/C0MCT/cV.T 

5ANPHS 


COMMON  HLOS ( 6  r 60 ) » RNG (6 . 60) .XANG(SOO) .YDB<300) . JP1 . JCNT » WlNDl » TOTL 
H»TOTLN*l51»lTST,NTSTrTHRML.THRD0,CD.ZBl,VFLAG.BANDL(ll) .FRE3<11) . 
206(11)  .BANDW.LOCaTO)  .DATE (2)  .AVELH30)  .  ASHI1  (30 ) . TDB ( 180) . TDEG  ( 18 
30) .KPl .TRECL.TLEV. TGDP.RTG.Q (300) » W (300  ) . TYTLE (8 ) » CAP1 (8) »CAP2( 13) 
4 ,C AP3 ( 13) »  CAP4 < 1 3 ) .CAP5Q3) »CAP61S) >CAP7<8) .NMS2.FREQ.PHINC 

DATA  RADcON/57. 2957795/ 

DATA  (SX(I)rI=lf4)/-1.0»-l. 0.1. 0.1,0/ 

DATA  (ST ( I ) .1=1. 4)/-l. 0.1.0. -1.0.1. 0/ 


COSPHS=CoS(PHIS) 

CV=CX/COSPHS 

CALL  VERTEX  (ZLO.ZHI) 

IF  (ZTG.LT.ZHI .OR.ZTG.GT.ZLO)  GO  TO  90 
TANPHS-SlN (PHIS) /COSPHS 
PHI=57.296*PHIS 

RAYTRACE  FROM  SOURCE  TO  BOTTOM 

CALL  RAYTRC  (PHIs.ZX ,ZBM,SXB.RXB,DRXB.PHlBOT.GRADB»ZEND) 
BRLrFUNU ( AHB » RADcON*PHIBoT . NUMH > 

RAYTRACE  FROM  SOURCE  UP  TO  SURFACE 

CALL  RAYTRC  (-PHiS.O.O.ZX.SXS.RXS.DRXS.PHlE.GRADS.ZEND) 
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ZSURFrO.O 

NS=0 

IF  (PHIE.NE.0.0)  GO  TO  10 
IF  (ZEND. NE. 0.0)  NS=1 
ZStJRF=ZEN0 
10  CONTINUE 

IF  (ZTG.EQ.ZX)  Go  TO  20 

RAYTRACE  FROM  SURFACE  (OR  VERTEX  DEPTH)  DOWN  TO  TARGET 
CALL  RAYTRC  (PHlE»ZSuRF»ZTG»SST»RSTrDRST.PHND*GRADT,ZEND) 

DIRS  IS  THE  NOISE  SOURCE  DIRECTIVITY 

EXPN=5.729*PHND 

IF  (PHND. LT. .349D  EXPN=2.0 

DIRS=1 .+(SIN(PHNd)**EXPN-1. ) * (2 . *FREG- . 02> 

IF  (ZEND.EQ.ZTG)  GO  TO  30 
WRITE  (6,100)  ZEND 
GO  TO  30 
20  SST=SXS 
RST=RXS 
DRST=DRXS 
PHND=PHIS 

30  CONTINUE 

R=2.0*RXB-RXS-R5t 

LIMIT  TO  FIVE  BOTTOM  BOUNCES 

IF.  (ICZ.EQ.0.AND.NMAX.GT.5)  NMAX=5 

DO  ,40  N— 1 #  NMAX 

EN=N 

HBOT=0.0 

IF  (ICZ.EQ.O)  HBoT=EN*BRL 
RN=2.0*EN*IRXB+RxS> 

SN=2.0*EN*(SXB+SXS> 

DN=2.0*EN*(DRXB+DRXS) 

DO  40  K=l»4 

PSs-SX (K) *PHIS 

IF  (PS.LT.PHMIN.oR.PS.GT.PHMAX)  go  TO  40 
RK=RN+SX(K)*RXS+sT(K)*RST 

maximum  range  check 

SK=SN+SX(K>*SXS+sT(K)*SST 
DK=DN+SX ( K ) *DRXS*ST ( k ) *DrST 

ELIMINATE  CAUSTIC  BY  RESTRICTING  OK 
IF  (ABS(DK).LE.l.)  Dk=1. 

HSB=HSURF(NfK)+HBOT 

CALL  HLOSS  (PHND, RK » SK * D« t GRADT r HK ) 

HK=HK+HSB+FUNU l ArESP , SX ( K  >  *PHI » NUMR ) 

CALL  SUMINT  (K  »N,NMAX»RK , HK  » PHI  *  DIRS) 

40  CONTINUE 


DIRECT  PATH 
41  K=5 
S=-1.0 

IF  (ZX-ZTG)  50»80.60 
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50  PN=PHIS 
PX=PHlS 
GO  TO  70 
60  PN=-PHIS 
PX=-PHlS 

70  RK3ABS<RXS+S*RST) 

if'(pn.lt.phmin.or.pn.gt.phmax)  goto  eo 

SK=ABS<SXS+S*SST) 

OK=ABS(DRXS+S*DRsT> 

HS=HSURFU*K) 

CALL  HLOsS  (PHND  ,  RK  r  SK »  Di\  »GRAOT  t HK ) 
HK=HK+HS+FUNU ( AReSP • PN»  NUMR ) 

CALL  SUMINT  (K»N,NMAX>RK,HKiPHI»DIRS) 
80  IF  (K.EQ.6)  GO  To  00 

SURFACE  REFLECTED  PATH 
K=6 

N=NMAX+2 

s=1.0 

PN=-PHIS 
PX=PHIS 
GO  TO  70 

90  CONTINUE 

UELR=R-RLAST 

RLAST=R 

return 

entry  preray 
NTST=0 
ITST=1 
RLAST=200.0 

nmax=bncs 

PHM AXzPHlD+O . 5+DeLPH 

phmin=phmax-oelph 

IF  (PHMIN.GE.  0.0)  ITST=2 
PHMAXzPHmAX/RADCoN 
PHMINrPHMlN/RADCoN 
RETURN 

100  FORMAT  (//5X,4HZeNDiE12.6//) 

END 


SUBROUTINE  RAYTRc  (PHlS»ZUPfZLO»SP»RPfDRDP*PHIEND»GRAD.ZNPl) 


PHIS  POSITIVE  MEANS  DOWN-GOING  RAY,  TRACES  FROM  ZUP  TO  ZLO 
PHIS  NEGATIVE  MEANS  UP-GOING  RAY,  TRACES  FROM  ZLO  To  ZUP 

COMMON  /COMA/  A VeLP ( 30 ) . AHB ( 30 ) *■ ARESP (30) »ASHIP(30> . ABW ( 30 ) , APROP ( 
130) , ADANF (30 ) ,NUMV,NUMH,NUMR,NUMs»NUMBW.NUmP.NUMF/COMCT/CV.TANPHS 


Rp=0.0 

SP=0.0 

DRDP=0.0 

SINPN=SIn(PHiS> 

COSPN=COs<PHlS) 

PHIM=PHIS 

N=1 
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C  TEST  FOR  NEGATIVE  PHlSTAkT 

IF  (PHIS. LT. 0.0)  GO  TO  20 

phis  positive*  Down-going  Ray 
zn=zup 

10  NUM=N 
LL  =  1 

GO  TO  40 

phis  negative*  up-going  ray 
20  ZN=ZLO 
30  NUM=NUMV-N 
LL=2 

40  AZN=AVELP(NUM) 

IF  ( (ZN.LT.AZN) .AND.  ILL. £0.1) >  GO  To  50 
IF  ( (ZN.GT.AZN) .AND, (LL.EQ.2) )  Go  TO  50 
N=N+2 

GO  TO  (10*30) *  Ll 
50  ZNP1=AZN 

IF  (ABS(ZNPl-ZN) ,LT. 0.001)  GO  TO  150 
IF  (ZNP1.GT.ZLO)  ZNP1=ZL0 
IF  (ZNP1.lT .ZUP )  ZNP1=ZUP 
CN=FUNU (AVELP  *ZN,NUMV ) 

CNP1=FUNU ( AVELP*  ZNP1 *  NUMV ) 

TEST  FOR  VERTEXInG  RAY 
IF  (CNP1.GT.CV)  GO  TO  60 
PHINP1=AC0S(CNP1/CV) 

IF  (PHIS. LT. 0.0)  PHINP1=-PHINP1 
C0SPN1=CNP1/CV 

SINPNI-SqRT ( 1 . 0-COSPNl*CoSPNl ) 

IF  (PHIS. LT. 0.0)  SINPNlr-SINPNl 


GO  TO  70 

VERTEXING  RAY  COMPUTATIONS 
60  ZNP1=ZN+(CV-CN)/(CNP1-CN)*(ZNP1-ZN) 
CNP1=CV 
PHINP1=0.0 
COSPN=1.0 
COSPN1=1.0 
SINPN1=0,0 
70  DELZ=ZNP1-ZN 

TEST  FOR  VERTICAL  PROFILE 
IF  (ABS(CNPl-CN) ,GT. 0.00001)  GO  TO  80 

C  VERTICAL  PROFILE  COMPUTATIONS 
DELRP=DELZ*COSPN/SINPN 
RP=RP+DELRP 
DELSP=DELRP/COSPn 
SP-SP+DELSP 
GO  TO  90 
C 

80  GRAD=(CNP1-CN)/DELZ 
rho=cv/grad 

DELRP=RH0*(SINPN-SINPN1) 
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KP=RP+DELRP 

DELSP=RH0* (PHIN-PH1NP1 ) 

IF  (DELSP.LT. 0.0)  ITEST=l 
SP=SP+DELSP 
C 

90  CONTINUE 

IF  (ITEST.NE.O)  WRITE  (6.160)  ZN.ZNPl *CN.CNPl»COsPNrCOSPNl .RP.SP 

itest=o 

TEST  FOR  LAST  TIME  THRU  FOR  VERTEXING  RAY 
IF  (ABS(CV-CNPI) ,LT. 0.00001)  GO  TO  120 

TEST  FOR  HORIZONTAL  ENTRANCE  RAY- 
IF  (ABS(PHIN) .LT.0.0001)  GO  TO  100 
DRDP=DRDP-DELRP'/SINPN/5INPN1 
GO  TO  110 

100  URDP=DRDP-RHo/ABs ( SINPN1 j 
110  CN=CNP1 
ZN=ZNP1 
PHIN^PHINPI 
N=N+2 

SINPN=SINPN1 
COSPN=COSPNl 

IF  ( (ZNP1.EQ.ZLO) .OR. (ZNP1.EQ.ZUP) )  GO  To  130 
IF  (PHIS)  30.10.10 

VERTEXING  RAY  -  MUST  EXIT 
120  CONTINUE 

DRDP=ORDP+RHo/ABs<SINPN) 

phiend^o.o 

GO  TO  140 

130  PHIENO=ABS(PHINPl) 

140  RPSRP/3000.0 
SP=SP/3000.0 
URDP=DRDP/3000.0*TANPHS 

return 

150  WRITE  (6,170)  ZNpl  ' 

RETURN 

160  FORMAT  ( 1X2HZN10X4HZNP18X2HCN10X4HCNP18X5HCOSPN7X6HCOSPN16X2HRP10X 
12HSP/8E12.6) 

170  FORMAT  (//*ZNP1=ZN=»E10.4) 

END 

SUBROUTINE  READIN  (WORD) 

THIS  SUBROUTINE  READS  THE  INPUT  PARAMETERS 

integer  pause.param.headeR.endaTa.prlos.word. BLANK. ASTRSK. DOT. VNAM 
1E.ANAME. SLSH.DAT A. EQU 

COMMON  /COMV/  VRbL ( 31 ) /COMA/ARR A Y ( 30 , 7) ,NUMA ( 7) /COMH/ I HDR ( 16) /COmX 
I/CX.ZBM.PHIC.ALPhAC.BION 

COMMON  HLOS (6.60) .RNG(6,60) ,XANG(300) » YDB ( 30 0 ) , JP1 , JCNT , WlNDI . TOTL 
ll.TOTLN. 151. ITST.NTST.THRML.THRDB.CD.ZBI, VFLAG, BANDL ( 1 1 ) ,FRE3(11>., 
2UB( 11) ,BANOW.LOCaT<3) . DATE ( 2 ) , AVEL1 (30 ) » ASHIl (30) .TDB(lBO) . TDEG( 18 
30) .KPi.TRECL.TLEv.TGDP.RTG, 0(300) r  W ( 300 ) ,TYTLE(8) »CAP1 (B) ,CAP2(13) 
4»CAP3(13) .CAP4(l3) »CAP5(13) fCAP6(8> ,CAP7(6) .NMS2.FREQ.PHINC 
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DIMENSION  VNAME ( 3 1 >  »  ANAmE ( 7 )  *  0ATA(5>,  lOATA(16)*  JDATA (81 ) 

DATA  (VNaME<I)*I  =  1s31)/5NZX  *5HZTg  *5HRMAX  *  5HPHID  * 5HDELPH , 5HB 

1NCS  *5H0UTPT,5HWSPD  , 5HDGREH , 5HSNPLT . 5HSSPLT t 5HALAT  * 5HBERNG  * 5HH0U 
2R  *5HACTiV*5HRAIn  . 5HHFLaG * 5hZTGl  * 5HZTGN0 * 5H2T INC * 5HRGT  .5HRGTNO 
3  * 5HRGI NC  »  5HT  ARG  * 5HT aRGO , 5HTaRG1 . 5HTARG2 * 5HTARG3 * 5HSHIPD .5HSHLFR * 
45HSHLFS/ 

DATA  NV/31/ 

DATA  <ANaME<I)*I=1*7)/5HAVELP*5HaHB  *5HaRESP*5HaSHIP,5HABW  r  5HAP 
1R0P  *  5H ADANF/ 

UAT A  NA/7/ 

UATA  PAUSE  » PAR AM, HEADER  * END AT  A  * PRLOS  * BLANK * ASTRSK t DOT/ ' PAUSE • * ' PA 
1RAM* * 'HEADE1 , *ENq-D' * 'PRLOS' ,  ' 

DATA  KOMA/' » '/»SLSH/'/'/,EQU/*= V 
DATA  KINN/1/.KOUT/2/ 

IF  (WORD. NE. PRLOS)  WRITE(6.230> 

READ  AND  PRINT  INPUT  CARD 
10  READ  (5*240)  K I Nn * ( IDAT A ( I X ) , IX=1 * 16 ) 

WRITE  (6*240)  KOuT * ( lOATA ( IX ) * IX=1 , 16 ) 

K=1 

DECODE (80*250. I  Data)  (JDaTA ( JX ) *JX=1*80> 

SEARCH  FOR  BEGINiNG  OF  LABEL 
20  DO  30  J=K*80 

IF  (JOATA( J)  .EQ.kOMA)  JDaTA(J)=BLANK 
IF  (JDATA(J).EQ.  SLSH  >  JDaTA.(  J)=BLANK 
IF  ( JDATA(J) ,NE. BLANK)  GO  TO  40 
30  CONTINUE 
GO  TO  10 

FIND  THE  END  OF  THE  LABEL 
40  DO  50  K=J*80 


if  (JDATA(K) .E&.EQU)  GO  TO  60 
IF  (JDATa(K) .EQ. BLANK)  GO  TO  60 
IF  (JDATA(K) .EQ.KOMA)  GO  TO  60 
IF  < JDATA(K) .EQ.SLSH)  GO  TO  60 
50  CONTINUE 
GO  TO  80 

PACK  UP  TO  FIVE  CHARACTERS  OF  LABEL  INTO  CONTROL  WORD 
60  K5=K 

JDATA(K)=BLANK 
IF  (K-J.GT.5)  K5=J+4 
DATA  (1)=BLANK 

ENCODE (10* 250*DATA)  ( JDATA ( JX) * DX=J * K5) 

WORD=OATA(l) 

IF  (WORD. EQ.PAUSt:. OR. WORD. EQ.ENDATA. OR. WORD. EQ. PRLOS)  RETURN 
IF  (WORD.EQ.HEADFR)  GO  To  210 
IF  (WORD.EQ.PARAm)  GO  to  20 
IF  (WORD. EQ. ASTRsK)  GO  To  200 

see  if  label  is  a  param  or  array 

DO  70  I=I*NV 

IF  (WORD.EQ.VNAMe( I ) )  GO  TO  90 
IF  (I.GT.NA)  GO  to  70 
IF  ( WORD ,EQ • ANAME ( I ) )  GO  TO  180 
70  CONTINUE 


64 


oonon  onoooo  no  o  o  no 


C  MUST  BE  AN  ERROR 
BO  WRITE  (6,260)  WOrD 
GO  TO  220 

SET  UP  TO  PROCESS  SINGLE  VALUED  PARAMETER 
90  ASSIGN  170  TO  JAd 
K1ND=1 

PUT  VALUE  string  together 
100  DO  110  J=K»80 

IF  (JDATA(J) .EQ. BLANK)  GO  to  120 
110  CONTINUE 
120  DO  130  K=J,80 

IF  (JDATA(K) .EQ.kOMA)  GO  TO  20 

IF  (JDATA(K) .EQ.sLSH)  GO  TO  20 

IF  (JDATa(K) .EQ.eQU)  JDATA(K)=BLaNK 
IF  (JDATa(K) ,NE. BLANK)  GO  To  140 
130  CONTINUE 
VAL=0.0 
GO  TO  ID 
140  NDOT=0 

DO  150  J=K.80 

IF  ( JDATA( J) .EQ«bLANK)  GO  TO  160 
IF  ( JDATA( J) .EQ.DOT)  NDOT-NDoT+1 
IF  (JDATA(J> .EQ.KOMA)  GO  TO  160 

IF  (JDATA(U) .EQ.SLSH)  GO  TO  160 

150  CONTINUE 
U=81 

160  KAR=JDATa(J> 

JDATA(J)=BLANK 

INSERT  A  DECIMAL  POINT  IF  NUMBER  HAD  NONE 
IF  (NDOT.EQ.O)  JdATA<J)=dOT 
DATA  <1)=BLANK 
DATA  (2)=BLANK 


ENCODE < 10 >250 .DATA) ( JDATA(JX) . JX=K, J) 

JDATA(J)=KAR 

FORM  THE  VALUE  As  A  REAL 
DECODE  (10.270, DATA)  VAL 
K=J 

GO  TO  STORE  PARAMETER  OR  ARRAY  ELEMENT 
GO  TO  JAD.  (170.190) 

STOR  SINGLE  VALUE  IN  THE  RIGHT  SPOT  AND  GO  BACK  FOR  MORE  IF  CARD  N 

170  VRBL ( I )=VAL 

IF  (K.GE.80)  GO  TO  10 
GO  TO  20 

SET  UP  To  PROCESS  ARRAY  ELEMENTS 
180  ASSIGN  190  TO  JAd 

VFLAG  IS  THE  VELOCITY  PROFILE  RERUN  FLAG  FOR  CURVaTuRE 

IF  (I.EQ.l)  VFLAg=0.G 

KIND=2 

NX=0 

GO  TO  100 
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C  UPDATE  number  of  elements  for  this  array  and  store  the  value 

190  NX=NX+1 

ARRAY(NX,I)=VAL 

NUMA<I>=NX 

IF  (K.GE.80)  GO  TO  10 
IF  (JDATA(K) .EQ. BLANK)  GO  TO  100 
GO  TO  20 

200  IF  (KIND.EQ.2)  Go  TO  100 
GO  TO  80 

READ  AND  STORE  THE  TITLE  INTO  THe  COMMON  AREA 
210  READ  ( 5 ,240 )  K INn r ( IHDR ( IX ) , IX=1 , 16 ) 

WRITE  (6,240)  KOuT » ( IHDR ( IX ) , IX=1 , 16 ) 

GO  TO  10 


220  CALL  EXIT 


230  FORMAT 
240  FORMAT 
250  FORMAT 
260  FORMAT 
270  FORMAT 
END 


( 1H1//) 


<T*.l6A5) 

C  S0A1 J 

(//lOX'ERROR  in  inputs 
(FIO.O) 


•  ,A5) 


SUBROUTINE  ScLNCe  (Tmin,tmax»nd»snmin»snmax, sninc) 

THIS  SUBROUTINE  PRODUCES  A  NICE  SET  OF  SCALE  NUMBERS  FOR  A  PLOT  AX 
INPUTS  ARE  THE  Data  min  »TMIN*»DaTA  max  *TmAX',  and  THE  REQUIRED  N 

intervals  along  the  axis  *nd‘.  the  routine  outputs  a  new  min  *snmi 
MAX  *SNMAX'»  and  increment  'SNINC*  FOR  NO+1  readable  values 

DIMENSION  POWTEN ( 7 ) *  TlCVAL(fi) 


DATA  (POwTEN(I)»I=l»7)/lo. ,100. YlOOO . » lOOOo . , 1 OOOOO , , lOODOOO . , 1000 

10000./ 

allowable  values  are  some  *ticval«  times  some  power  of  ten 
DATA  (TICVAL(I)»I=1»8)/ 0.5,1. 0, 1.5,2. 0,5. 0,10. ,15. ,20./ 

uiv=nd 

set  up  rounding  value 
ROUNOrO.ol 
st=tmin 

IF  (ST.LT.0.0)  RoUNDsROUND-I.O 

DETERMINE  first  GUESS  at  increment 
0=(TMAX-TMIN)/DIv 
P=1.0 

IF  (D.LT.10.0)  Go  To  30 

c  reouce  increment  to  range  d  -  lo) 

DO  10  1=2,7 

IF  (POWTEN(I) .gt.d>  GO  To  20 
10  CONTINUE 
20  P=P0WTEN(I-1) 

30  D=0/p-0.01 
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C  FIND  FIRST  ALLOWABLE  TICVAL 

00  40  1-1 *7 
IS=I 

IF  (TICVAL(I) .GE.D)  GO  To  50 
40  CONTINUE 

compute  new  INCREMENT 
50  D=TICVAL< IS)*P 

ST=D*AINT(ST/D+RoUND) 

TEST  NEW  increment  To  make  sure  data  will  fit  in  new  range 
TEST=ST+<D1V+0.0i)*O 
if  (TEST.GE.TMAX)  go  to  60 

is=is+i 

GO  TO  50 

compute  new  min  and  adjust 
60  ST=ST-AlNT< IDIV+(ST-TMAX)/D)/2.0)*D 
IF  (ST*TmIN>  70*70*  80 
70  ST=0.0 
80  SNMIN-ST 

COMPUTE  NEW  MAX 

snmax=st+d*div 

SMNC=D 
RETURN 
end 


SUBROUTINE  SHIPlN  (SHIPI) 

THIS  SUBROUTINE  CONSTRUCTS  THE  SHIPPING  DENSITY  NOISE  INTENSITY 

histogram 


COMMON  /COMV/  ZX , ZTG. RMAX  *  PHID*  DELPH* BNCS  * OUTPT  *  WSPD  * DGREH » SNPLT *  S 
ISPLT  *  ALAT  *  BERNG  *HOUR *  ACT  IV  * RAIN*  HFLAG  *  ZTGl » ZTGNO  *  ZTINC  *RGT  *  RGTNO  *  R 
2GlNC*TARG*TARGO*TARGl*TARG2*TARG3.SHIPD*SHLFR,SHLFS/COMA/AVELP(30) 
3*AHB(30) ,ARESP(3o)*ASHIP130)*ABW(30) * APROP(30)*ADANF(30) *NUMV*NUmH 
4 , NUMR  *  NUmS  *  NUMBW , NUMP  *  NUmF 

COMMON  HLOS(6*60) *RNG (6*60) *XANG( 300) *YDB (300) »JP1*JCNT*WINDI*T0TL 
HfToTLN* 151  * ITST, NTST*THRML» THROB* CD rZBi* VFLAG*BANDL ( 11 ) *FRE3(11)  * 
20B(11)_*BaNDW,L0CaT(3) *  DATE (2) ♦ AVELI (30) * ASHI 1 ( 30 ) ,TDB(1B0) *  TDEG ( 18 
30) *KP1 *TreCL*TLEV*TGDP*RTG*Q(300) *W<300) .TYTLE(S) *CAP1 (8) »CAP2( 13) 
4  f CAPS ( 13 ) *CAP4 ( 13) *CAP5(i3) *CAP6(8) *CAP7(8) *NMS2*FREQ*PHINC 

REAL  LMULT 


J=0 

K=NUMS-2 

DO  40  I=1*K*2 

INITIALIZE  PARAMETERS 

LMULTrl.o 

SMULT=1*0 
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n  n  n  r>  r> 


c 

c 


10 

c 

c 


20 


C 

C 


30 


c 

c 


STORE  SHIPS  AND  FIRST  RANGE  OF  RANGE  INCREMENT 

RANG1=ASHIP(I> 

SHIPS=ASHIP<I-H> 

J=J+2 

IF(ASHIP< 1+2) .GE.0.0)  GO  TO  30 

CALCULATE  INTENSITY  MULTIPLICATION  FACTORS  FOR  SHIP  SPEED 
AND  LENGTH 

lF(ASHIP(I+2) .LE.-50.)  GO  To  20 
ADSHP=50.*ALoG10(-ASHIP<1+2)/12.) 

SMULT=2 . ** ( ADSHP/3 . 01 03 ) 

1  =  1  +  1 
GO  TO  10 

ADLNG=2O.*ALOGl0(-ASHIP<I+2)/300.) 

LMULT=2.**UDLNG/3.0103) 

1  =  1+1 
GO  TO  10 

calculate  area  And  construct  shipping  density  noise  intensity 

HISTOGRAM 

RANG2=ASHIP( 1+2) 

AREA= . 0087266* (RANG2**2-RANG1**2) *DGREH 

dens=ships*ioooo./area 
ash i i ( j_i ) -Rang 1*2 « 

ASH II (U)=DENS*SHiPI*SMULT*LMULT 

continue 

NMS2=J+1 

ASHI1 (J+1)=RANG2*2. 

RETURN 

END 


subroutine  SSPLOT 

THIS  SUBROUTINE  wILL  PRODUCE  A  CALCOMP  PLOT  OF  A  SOUND  SPEED  PROFI 

•disspla*.  the  Depth  values  are  along  the  x  axis  and  the  vel  value 
the  y  axis,  the  plot  will  be  DRAWN  on  a  a.5  X  11  INCH  page  with  au 

DIMENSION  LABEL (l4) »  X(15)»  Y(15) 

common  /com v/  zx,ztg.rmax»phid»delph. bncs » outpt . wspd. dgreh. snplt ,s 

1SPLT . ALAT  » BERNG  *  HOUR. ACT I V  »RAIN>  HFLAG » ZTGl » ZTGNO. ZTlNC . RGT*  RGTNO .R 
2gINC  f TARG  t TARGO . TARG1 *  TarG2 i T ARG3 . SHIPD . SHLFR . SHLFS/COMA/A VELP <  30 ) 

3. AHB<30)  .AREsPOo)  »ASHIP(30)  ,ABW(30)  .  APROP  <  30  > » ADANF  ( 30 ) » NUMV ,  NUMH 

4 ,  NUMR » NUMS . NUMBW , NUMP f NUmF 

COMMON  HLOS ( 6 • 60 ) »RNg(6»60) ,XANG<300) .YDB<300) . JP1 . JCNT.WINDI »TOTL 
H.T0TLN.I51. ITST,NTST.THRML.THRDB»CD.ZB1.VFLAG.BANDL(11) .FRE3(ll) , 
2DB (11) r BANDW. LOC aT (3) » DATE (2) i AVeLI < 30 ) r ASHI 1 ( 30) »TDB<180> .TDEG<lfl 
30) .KPlrTRECL.TLEv*TGDP»RTG.Ql300) ,W(300).TYTLE(B) .CAP1 ( 6 ) »CAP2 C 13) 
4  r  CAP3 ( 13) >  CAP4 ( 13) .CAP5 ( l3) .CAP6(8) rCAP7<8) »NMS2,FREQ»PHINC 
DATA  NP/100/ 

C 

NP=NP+1 

NV=NUMV/2 
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C  BEGIN  PLOT 

CALL  BGNPL  (NP) 

IF  (HFLAG.LT.O.O.OR«SNPLT.GT.O.O)  GO  TO  10 
READ  ( 5  r  50  J  (LOCaT(I) *1  =  1*3) 

READ  (5.50>  (DATe(I) ,1=1,2) 

10  CONTINUE 

ENCODE  (40,  7o»CAp2)  (LOCATd)  .1  =  1,3) 

ENCODE (40,  80  » C Ap3 ) (DATE(I> , 1  =  1 '2) 

CALL  TITLE  (1H  ,-l , 0 , 0 . 0 . 0 , 5 .0 * 3 . 0 ) 

C 

C  FIND  RANGE  OF  SOUND  speed 

YMIN=6000.0 
YMAX=0  •  0 
DO  20  1=1, NV 
X(I)=AVELK2*I-1) 

Y(I)=AVEL1(2*I) 

IF  (Y(I) .GT.YMAX)  YMAX=Y(I) 

IF  <Y<I) .LT.YMIN)  YMlN=Y(I) 

20  CONTINUE 
C 

C  FIGURE  OUT  NICE  SCALES  FoR  X  AND  Y 

CALL  SCLNCE  (YMIn»YMaX,5. YMN,YMX,YINC) 

CALL  SCLNCE  ( X (1 ) r X ( NV ) , lO , XmIN , XMAX , X INC ) 

C 

C  SET  UP  AND  DRAW  THE  X  AND  Y  AXIES 

CALL  BANGLE  (160.0) 

CALL  BSHIFT  (4. 5, 3.0) 

CALL  GRAF  (XMIN,XINC,XMAX,YMN,YINC,YMX) 

C 

c  set  up  the  required  alphabet 

CALL  BASALF  ( 'L/cSTD* ) 

CALL  MIXaLF  ( 'STANDARD’ > 

CALL  HEIGHT  (0.1) 

C 

CALL  CURVE  (X.Y.nV.O) 

C 

c  make  and  label  X  AXIS 
xtmp=xmin 

DO  30  1=1,11 
IXTM=XTMp 

ENCODE  (10,60*  LAbEL ( 13) >  IXTm 
LABEL(I)=LABEL(l4) 

30  xtvp=xtmp+xinc 
CALL  XAXANG  (90.) 

CALL  XLBAXS  (LABEL, 1, 11.5.0, »  $*,100,0.0,0.0) 

C 

C  MAKE  AND  LABEL  Y  AXIS 
YTMP=YMN 
DO  40  1=1,6 

xytm=ytmp 

ENCODE  (10, 60, LAgEL (13) >  IYTm 
LABEL ( I ) zLABEL ( I4 ) 

40  YTMP=YTMP+YINC 

CALL  YLBAXS  (LABEL , 1 , 6 , 3 . 0 , • (SOUND  SPEED)  -  FT/SEC$ * , 1 00 , 0 . 0 » 0. 0 ) 
CALL  ANGLE  (-160.) 

CALL  MESSAG  (MDePTH)  -  FEET ' , 14 , 3 , 0 ,-0 . 6) 

CALL  ANGLE  (90. ) 

CALL  RESET  ( 'BASaLF' ) 

CALL  MESSAG  (C AP2 , 100 * -1 .25 , .45 ) 

CALL  MESSAG  ( CAP3* 100 , -1 . 0 , . 45 ) 
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FRAME  AND  END  PLOT 
CALL  FRAME 
CALL  ENDPL  (0) 
RETURN 


50  FORMAT  (3A5) 

60  FORMAT  (110) 

70  FORMAT  C »LOCATIOn‘ »4Xr3A5# 'S' ) 
80  FORMAT  ( *DATE*#l3X»2A5. ) 

END 


SUBROUTINE  SUMINT  (KiNrNMAX»RK»HK»PHI pDIRS) 

THIS  SUBROUTINE  COMPUTES  THE  NOISE  INTENSITY  ALONG  EACH  RAY 

COMMON  /COMV/  2 X , ZTG . RM Ax  *  PHI D t DELPH , BNCS » OUTPT t WSPD » DGREH r SNPLT , S 
lSPLT,ALAT.BERNG»H0UR,ACTlV,RAlN.HFLAG»2TGlfZTGN0,ZTINC»RGT»RGTN0tR 
2GlNC-TARG»TARGO»TARGl»TARG2*TARG3»SHIPDrSHLFR»SHLFS/COMA/AVELP(30) 
3f AHB<30>  » ARESP(3o) , ASHIP<30) , ABW ( 30 ) . APROP ( 30 ) ,ADANF(3Q) .NUMVfNUMH 
R.NUMR,NUmS»NUMBW.NUMP*NUmF/CQMX/CX.ZBM»PHIC* alphac*bion 
COMMON  HLOS(6»60) »RNG<6»60) *XANG(300) »YDB(300) » JP1 » JCNT»WINDI »TOTL 
HfT0TLN#I51rI TST  , NTST  *  THRML  > THRDB  f CD . ZB1 »  VFLAG  *  BANDL ( 11 ) ,FRE3(11) > 
2DB (11) »BANDW»L0CaT(3) t DATE (2) »AVEL1(30) fASHll(30) fTDB(180) »TDEG(18 
30) »KPl#TRECL,TLEV»TGDP»RTGrQl300) »W(300) »TYTLE<8) ,CAP1 < 8) *C AP2 < 13 ) 
4  » C AP3 (13) rCAP4(l3) rCAP5(l3) »CAP6(8) »CAP7(8) »NMS2,FREQ»PHINC 

DATA  RADC0N/57. 2957795/ 

TEST  FOR  FIRST  BOUNCE 

IF  (NTST.EG.l)  Go  TO  10 

NTST=1 

DOWNI=0.0 


UPI  =  0.0 

10  ABSFI=ABS(PHINC) 

RG=RK 

SHIP1 =FUNS ( ASH  1 1 , RG »  NMS2 ) 

DR=ABS(RNG (K»N)-RK) 

IF  (RK.Eq.O.O.OR.RNG(K*N) .EQ.0.0)  DR=5,0 
HL=.5*(HK+HL0S(K,N) ) 

HLOS ( K  f  N ) =HK 
RNG(KiN)=RK 

TEST  FOR  FIRST  T iME  THRU 
IF  (ITST.EQ.I.ANd.K.EQ.6)  GO  TO  70 

IF  (ITST.EQ.  2  -AND.  K  .EO.  2  .AND.  N  •  EQ  •  NMAX)  GO  TO  70 
IF  (ITST.GE.l)  Go  TO  50 
AREA=DGREH*.00873*DR*«2.0*RG+DR) 

XZZ  =  ARE  A»0.5**(H|_/3. 0103) 

ASSUME  WIND  NOISe  DIPOLE  DIRECTIONALITY  PATTERN  AT  HIGH  FREQ 
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nnoo  r>  n  on  n  n  o 


20 


30 


1*RECJ=X2Z*WI.NDI*DIRS 
SRECI=XZZ*SHIP1 
T0TI=WRECI+SRECI 
IF  (K.GE.  3)  GO  TO 
DOWNI=DOWNI+TOTI 
GO  TO  30 
UPI=UPI+TOTI 
PHUP=-PHI 
CONTINUE 

IF  (K  .Eq.  2  .AND. 


20 


I TST  .EQ.  -1  .AND.  N  .Eg.  NMAX)  GO  TO  31 


IF  (K.NE.6)  GO  To  50 
31  CONTINUE 
NTST=0 

OOWNI=DOwNI+(THRmL+BION)*ABSFI 
UPI=UPI+(THRML+BiON)*ABSFI 
TOTLI=TOTLI+DOWNi+UPI 
DBDN=10.Q*ALOG10 (DOWNI )+ 184.8856 
DBUP=10.0*ALOG10{UPI >+184.8856 

NORMALIZE  LEVEL  TO  ONE  STERAD  AND  STORE  ARRAYS  FOR  NOISE  PLOT 


JCNT=JCNT+1 

CALCULATE  PER  STeRADIAN  CORRECTION  FACTOR 
AODB=-10.0*ALOG1o<ABSFI*cOS(PHI/RADCON >4.1 09662) 

J=2*UCNT-1 

JP1=J+1 

XANG(J)=PHUP 

YDB(U)=DBUP+ADDB 

XANG( JP1>=PHI 

YDB ( JP1 ) =DBDN+ADdB 

40  CONTINUE  „ 

IF  (SNPLT.GT.0.0  .OR.  OUTPT  .LE»  0.0)  GO  To  50 

PRINT  90,  PHl»PHlNC*-RNG<4»NMAX) »HL0S(4,NmAX> »OBUP#DBDN#ADDB 
50  RETURN 


70  ITST=0 

IF  (K  .Eq,  2)  ITST=-1 
RETURN 


90  FORMAT  (16X»8F12.3) 

END 

SUBROUTINE  SURFI 

THIS  SUBROUTINE  CALCULATES  THE  INTENSITIES  OF  ALL  THE  SURFACE 
NOISE  GENERATORS 

COMMON  /COM V/  ZX , ZTG , RMAX rPHID> DELPH, RNCS * OUTPT, wSPD , DGREH* SNPLT, S 
lSPLTiALAT*BERNGrHOUR, ACTIV,RAIN*HFLAG» ZTGl*ZTGNO,ZTlNC.R&T,RGTNO,R 
2GlNC,TARGiTARG0»TARGlpTARG2*TARG3*SHIPD,SHLFRf SHLFS/COMA/AVELP(30) 

3 ,  AHB ( 30 ) ,ARESP<30>  *ASHIP(30) »ABW<30) ,APROP(30) , ADANF (30 ) »  NUMV t NUMH 

4 ,  NUMR , NUmS • NUMBW , NUMP  r  NUmF/COMX/CX , ZBM , PHI C • ALPH AC , B ION 

COMMON  HLOS (6  f  60 ) , RNG C 6 , 60 ) r XANG ( 300 ) r YDB ( 300 ) , JP1 , JCNT » to INOI t TOTL 
H»T0TLN»l51»lTST,NTST»THRML,THRDB*CD,ZBl,VFLAGrBANDL(ll) ,FRE3(11) , 
2DB(11) *BANDW,LOCAT(3) ,DATE(2) »AVeL1(30),aShI1(30) ,TDB(180) ,TDEG(18 
30) , KPl r  TrECL  >  TLEv  *  TGDP  »RTG»Q(3Q0) , W ( 300 ) ,TytLE(8) ,CAP1 (8) ,CAP2 ( 13) 
4 ,C AP3 ( 13 ) ,C AP4 (13 ) »CAP5<l3) ,CAP6(8) ,CAP7(8) , NMS2 , FREQ , PHINC 
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WlNT*  RAINI  AND  aSHIP(2N)  ARE  SOUND  INTENSITY  DENSITIES  DUE 
TO  WIND*  RAIN  ANd  SHIPS.  •  {WATTS/M**2)/KYDS**2 

DO  10  J=l*ll 
BB1=AL0S10(FRE3(j) ) 

10  DB(J)=B81* (4.22*bB1-33.4J-1 0.9 
CALL  BWIDTH  (THRmL) 

ESTABLISH  LOWEST  THERMAL  LEVeL/VeRT  DEG/HOR  DEG 

THRML=THRML*DGREH*1 .5432E-5 
THRD6=10.0*ALOGl0(THRML+BlON)+84.8956 

ESTABLISH  MINIMUM  5CALE  CD  FOR  NolSE  PLOT 
AND  CONVERT  TO  Db/MICROPaSCALS 

CD=?0. 

IF  (THRDb,LE.-B0.0>  CD=10. 

IF  (THRDB.LE.-90.G>  CD=0.0 
IF  (THRDB.LE.-lOo.O)  CD=-10,0 
IF  (THRDB.LE.-11o.O>  CD=-20.0 

WlNT=0.0 

IF  (WSPD.EQ.0.0)  GO  TO  30 

WMULT  IS  WIND  SPeED  DEPENDENCE 

WMULT=.065*WSPD**.7 
DO  20  J=1*U 
BB1=ALOG10(FRE3(J) ) 

CORW  IS  THE  SLOPe  and  LEVEL  CORRECTION  For  SURFACE  NOISE  SPECTRA 
EXCEPT  SHIPPING  NOISE 

C0RW=29.6*10.** Ii.43E-5*fRE3(J>**1.54) 
wLEVl=-la.*BBl-4, 

WLEV2=9.66*BB1*Bb1-43.99*BB1+23.33 
20  DB(J)=WLEV1+(WLEV2-WLEV1)*WMULT+C0RW 


CALL  BWIDTH  (WIND 
30  RAINI=0.0 

IF  (RAIN, EO. 0.0)  GO  TO  50 

CALCULATE  NOI5E  DUE  TO  RAIN 

ARAINrALOGlO (RAIN) 

DO  40  J=l*ll 
BB1=ALOG10(FRE3(J) ) 

COR W=.0667*BB1**4. 106+29. 64 
40  DB(J)=C0RW+5.5*Bb1+14.5*aRAIN-49.5 
CALL  BWIDTH  (RAInU 
50  wINDI=WINT+RAINI 

CORL  IS  THE  MODEL  DISTORTION  FACTOR  FOR  OMNI  SOURCES 

IF  (SHIPD.EQ.0.0)  GO  TO  90 
DO  60  J=l*ll 

BBl=ALOG10(FRE3(j)-2.*SHlPD+12. ) 

C0RL=14.*FREO+8.74+4.*SHlPD 

SL1=-1.0*BB1+I.l6 
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SL2=3.3*BBl-6.274 

DB< J)=-l5»-10.*ALOG10 ( 1 0 . **SLl+l 0 . **SL2 ) +CORL 
IF(SHLFS  .LE.  0.0)  D6( J)=DB( J)+5. 

60  CONTINUE 

CALL  BWIdTH  (SHIpI ) 

70  CONTINUE 

construct  shipping  density  NOISE  INTENSITY  HISTOGRAM 
CALL  SHIPIN  (SHIPI) 

RETURN 

90  SHIPI=0.0 
GO  TO  70 
END 


subroutine  target 

THIS  SUBROUTINE  COMPUTES  THE  DIRECTIONAL  SIGNAL 
LEVEL  RECEIVED  FrOM  A  TARGET 


COMMON  / C0MV/  ZX,ZTG,RMAX*PHIDr DELPH,BNCS'OUTPT,wSPD.DGREH.SNPLT.S 
lSPLT»ALATrBERNG>HOUR.ACTiV.RAIN*HFLAG>ZTGl.ZTGNO»ZTlNC.PGT»RGTNO,R 
2GINC»TARG»TARG0*TARG1»  TARG2  »  TARG3  *  shipo. shlfr  »  SHLFS/COMA/AVELP ( 30) 

3  t AHB ( 30 ) ,ARESP(3o) , ASH IP (30) »  ABW ( 30 } rAPRoP<30) »  ADANF ( 30 ) *  NUMV  »  NUmH 

4  f NUMR»  NUMS»  NuMBW , NUMP  t NUMF/COMX/CX » ZBM» PHlC  » ALPHAC  t BION/CQMCT/CV » T 
5ANPHS 

COMMON  HLOS (6*60) ,RNG(6.b0) , XANG1300) .YDB<300) »  JP1 »  JCNT  #  W2  NO  I ^TOTL 
II pTOTLN* I5l> ITST,NTSTf THRMLr THROB* CD fZBl. VFLAG ► BANDL < 1 1 > fFRE3(U) t 
2DB< 11) .BANDW»L0CaT(3) . DATE ( 2) » AVeLI ( 30 ) » ASHI 1 ( 30 ) , TdB ( 180 ) rTDEG<18 
30) »KPl.TRECL»TLEV'TGDP»RTG»Q(300) »W<300) ,TYTLE<8) ,CAP1<6) *CAP2(13) 
4.CAP3(13) ,CAP4(13) » CAP5 ( 13 ) . CAP6 ( S ) . CAP7 ( 8 ) ►NMS2,FREQ .PH  INC 

DATA  RADC0N/57, 29577/ 


CX=FUNU(AVELP.ZX,NUMV) 

CTG=FUNU(AVELP*TgDP»NUMVJ 

PHITG=0.0001 

IF  (CX.LE.CTG)  PhITG=ACOS(CX/CTG) 
TDEG(l)=90. 

TDB(l)=-300. 

KP1  =  1 

UP=1.0 

RUN=0.0 

YLASTrO.O 

TRECI=1.0E-30 

PRINT  230 

IF  (TARG.LT. 0.0)  GO  TO  20 

calculate  broadband  target  energy 


DO  10  J=l»ll 
BB1=ALOG10 (FRE3(j> ) 

10  QB( J)=TARG0+TARGl*BBl+TAHG2*bBl**2+TARG3*BBl**3 
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CALL  BWIDTH  (TRGtD 
GO  TO  30 

20  TRGTI=3.2467£-9*i0.**(TARG0/10.  > 

30  CONTINUE 

TLEV  IS  EFFECTIVE  TARGET  LEVEL  IN  DB 

TLEV=10.*ALOG10  <TRGTl)+84„S8b6 

FIND  BOTTOM  GRAZING  ANGLE 
PHIBH=ACOS(CX/AVELP(NUMV) ) 

40  PHI1=DELPH/(2.0*RADCQN) 

IF  (RGT.GE.500. )  GO  TO  200 
CZ=-1.0 
ITEST=0 
UPHI=-.05 
TINTl=1.0E-35 
50  IR=0 

XTRAB=0.(j 

TEST  FOR  RAY  WIThIN  ACCEPTABLE  ANGULAR  LIMITS 

IF  (PHII.GT.PHIBm.AND.CZ.GT.0,0)  GO  TO  lao 

IF  (PHI1 .LT. 0.0. OR. PHI1.lT.PHIBM.AND.CZ.lT. 0.0)  GO  TO  200 

COSPl=COS(PHll) 

CV=CX/COSPl 

CALL  VERTEX  (ZLO.ZHI)  , 

TEST  FOR  RAY  UNAbLE  TO  REACH  TARGET  DEPTH 

IF  (ZLO.LT. (TGDP  +  5. ) .OR.ZHI ,GT. (TGDP-5,  ))  GO  TO  160 
TANPHS=TAN(PHI1> 

CALL  RAYTRC  (-PHi 1 , 0 . 0 r Zx » SI  * R1 » DRDP1 . PHnDI » GRADl » ZEN01 5 
CALL  RAYTRC  (PHNol »  ZENOi  ,  TGDP  .  S2  .  R2  .  DRDP2  •  PHND2  .  GPAD2  »  ZEIMD2  ) 
CALL  RAYTRC  (PHI  1 . ZX ,ZBM. S3* R3 . DRDP3 » PHND3 , GRAD3 . ZEND3 ) 

TEST  FOR  RELATIVE  POSITION  OF  TARGET 

IF  (TGDP-ZX5  60» 60.70 


TARGET  ABOVE  RECEIVER 
60  RG1=2.*R2 

RG2=2  ,* (R3+R1-R2) 

RSUP=Rl-R2 

RSDN=Rl-R2+2.*R3 

Elsi.O 

E2=0.0 

GO  TO  80 

target  below  receiver 
70  RG1=2.*(R3+Rl-R2) 

RG2=2,*R2 
RSUP=Rl+R2 
RS0N=R2-R1 
£1  =  0.0 
£2=1.0 

TEST  FOR  UP-GOINg  OR  DOWN-GOING  RAY 
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80  IF  (UP)  90  »  90  r 10o 
C  DOWN  GOING  RAY 

90  R5TRT=R50N 
Fl=1.0 
F2=0.0 
GO  TO  110 

C  UP  GOING  RAY 

100  RSTRT=RSUP 
F1=0.0 
F2=1.0 

110  5IZE-RG1+RG2 

CALCULATE  number  OF  CYCLES  TO  target 

IF  (ITEST.EO.O)  N=<RTG-RSTRT)/SIZE 
Y0-N 

rnode=yo*size+rstrt 

ITESTrl 

TEST  TO  GET  THE  SAME  NODE  AS  FIRST  TIME  THRU  AND  LOCATION  OF 
TARGET  IN  CYCLE 

IF  (IR)  130  « 120  . 140 
120  IF  (RNODE+RGI-RTg)  130.130.140 
130  RN0D1=RN0DE+RG1 
Gl-1  •  0 
IR=-i 
GO  TO  15Q 
140  RNOD1=RNODE 
G1=0.0 
IR=1 

150  CONTINUE 

TEST  FOR  RAY  WITHIN  RANGE  INCREMENT 

DR22RN0D1-RTG 
IF  (DR2)  160.170,170 

160  PHI1=PHI1+DPHI*<1.0+PHI1*.13*C1.+CZ> ) 

GO  TO  50 

CALCULATE  LOSSES 
170  CONTINUE 
ITEST=0 

SNODl=2.*Y0*(Sl+s3)+El*<Sl-S2+Fl*2.*S3>+E2*(S2+Sl*<F2-Fl) )+2.*Gl*( 
1E1*S2+E2*(S3+51-S2>  ) 

FNOD1=2.*YO*(DRDP1+DRDP3)+E1*(DRDP1-DRDP2+F1*2.*DRDP3)+E2* (DRDP2+0 
1R0P1* (F2-F1 ) )+2.*G1*(E1*uRDP2+E2*(DRDP3+DRDP1-DRDP2) ) 

CALL  HLOSS  ( PHND2 .RNoDl »  SNOD1 » FNQD1 . GR AD2  »  HLOS1 ) 

TEST  FOR  BOTTOM  BOUNCE 

BLOS=O.Q 
SLOS=0.0 

IF  (ZEND3.GE.ZBM)  BLOS=FuNU < AHB . PHND3*R ADCON » NUMH) 

IF  (ZEND1.LE.0.O)  SLOS=YO*.5 

IF  (UP.LE.0.0.AND.TGDP.LE.ZX.OR.TGDP.GE.ZX.AND.G1.EQ. 1.0)  XTRAB=1, 
10 

PHI=-uP*PHIl*RADcON 

TLOS= ( Y0+XTR AB) *BLOS+HLOS1+FUNU ( ARESP .PHI » NUMR ) 
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ooo  o  ooo  non  non 


sum  intensity  OVER  one  degree  and  store 

KPl=KPl+l 

TDB ( KP1 ) zTLEV— TLqS-SLOS 
TDEG(KP1)=PHI 
TGINT=10.**<TDB(kP1)/10. ) 

TINT1=TInT1+TGINj 

TRECI=TReCI+TGINT 

IF  (ABS<TDEG(KPl)-TDEGtKPl-l) ) .Lt.1.0)  GO  TO  210 
IF  ( YO.NE. YLAST. AND. Y0 .GT.5. )  TInT1=2. *TlNTl 
TDB(KP1)=10.*ALOg10 (TINT!) 

PRINT  24a#  TDEG(kPI) .TDB(KPl) 

TlNTl=1.0E-35 

IF  (RNODI.EQ.RSTrT)  go  TO  180 

TEST  FOR  CHANGE  FROM  UP  TO  DOWN  GOING  RAYS 

IF  (PHII.GT.PHIBm- AND.CZ. LT. 0.0. OR. PHI l.LT. PHI BM. AND.C2. GT. 0.0)  GO 
1  TO  110 

IF  (CZ.GT.0.0)  Go  TO  180 
GO  TO  200 
180  RUN=RUN+1.0 
UP=-1 .0 
ITEST=0 

IF  (RUN-1.0)  40  »  40  » 190 
190  TRECL=10.*AL0G10(TRECI) 

PRINT  220#  TRECL,ZX#TLEV,FREQ»TGDP#RTG 
RETURN 

START  bounce  free  rays 

200  PHI1=PHITG 
CZ=1.0 

DPHI=2.5/RTG 

IF  (OPHI.GT..005)  DPHI=.005 

itest=o 

TINTl=1.0E-35 
GO  TO  50 
210  KPl=KPl-l 
YLA5T=YQ 
60  TO  110 

220  FORMAT  (//F12.1#»  DB  SIfi  REC  AT'.FT.O#'  FT  FROM*  »F7.1#*  DB  TARGET 

1  AT • # F6 • 4 #  *  KHZ #  *  # F6. 0 »  *  FT  DEPTH# * #F8.3# *^KTD*  /lMl) 

230  FORMAT  (/30X.*REC  ANGLE 5X# *TARGET  LEVEL  * /33X .* DEG '# 10X #♦ DB/M  PA 
W) 

240  FORMAT  (2BX » F10.3 » 5X . Flo  .2 ) 

END 


subroutine  vertex  (zlo.zhi) 

THIS  SUBROUTINE  COMPUTES  HIGHEST  AND  LOWEST  DEPTHS  REACHED  BY  RAY 

COMMON  /COMV/  ZX,ZTG.RMAX#PHlD#DELPHrBNCS#OUTPT# WSPD#DGREH.SNPLT#S 
1SPLT#ALAT#BERNG#H0UR,ACTIV#RAIN#HFLAG#ZTG1#ZTGN0#ZTINC»RGT#RGTN0,R 
2GlNC#TARG»TARGO#TARGl#TARG2fTARG3»SHIPD#sHLFR#SHLFS/COMA/AVELP(3o) 
3»AHB(30) ,ARE5P(3o) #ASHIP(30) #ABW(30) #APRoP<30>  rADANF (30) #NUMV#NUMH 
4  »  NUMR #  NUMS  #  NUMBW , NUMP  #  NUmF/COMCT/C V . TANPHS 
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o  o 


c 

N=-l 

10  N=N+2 

ZN=AVELP(N> 

IF  (2N.LT.ZX)  GO  TO  10 
NST=N-2 

SEARCH  UP 
ZNP1=AVELP(N) 
CNP1=AVELP(N+1) 

N=NST 

20  ZN=AVELP(N) 

CN=AVELP(N+1) 

IF  (CN.GE.CV)  GO  TO  30 
CNP1=CN 
2NP1=ZN 
N=N-2 

IF  (N.GT.O)  GO  To  20 
ZHI=0.0 
GO  TO  50 

30  GRAD=(CN-CNP1)/<ZN-ZnP1> 
IF  (GRAD. EQ. 0.0)  GO  TO  4u 
ZHI=ZNPl+(CV-CNPl)/GRAD 
GO  TO  50 
40  ZHI=ZNP1 
C 

C  SEARCH  DOWN 

50  CONTINUE 

ZN=AVELP(NST) 

CN=AVELP(NST+1) 

N=NST+2 

60  2NP1=AVELP(N) 
CNP1=AVELP(N+1) 

IF  (CNP1.GE.CV)  GO  TO  70 

CN=CNP1 

ZN=ZNPl 

N-N+2 

IF  (N.LT.NUMV)  Go  TO  60 
-  ZL0=AVELP(N-2) 

return 

70  GRAD=<CNP1-CN)/<ZNP1-ZN) 
IF  (GRAD.EQ-O.O)  GO  TO  flu 


ZLO=ZN+(CV-CN)/GRAD 

return 

80  ZL0=ZN 
RETURN 
C 

END 
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